Read in individual subject rating data.
Columns are:
id (subject ID, possibly recode for public release)
Age (in years)
Sex (Male or Female)
InjectionUser (TRUE for participants who report a history of injectin methamphetamine, opioids, or both)
ImageSet: control, opioid, or meth
File: file name
Category: subcategory within ImageSet. Possible categories are:
ImageSetFile:
Order: the order the images were seen for this subject, from 0 to 359, spans across two visits.
Visit: V1 or V2 for the two visits
time: the order the images were seen for an individual visit, from 0 to 179
RatingType (which question this response is for: valence, arousal, craving, related, typicality)
Questions were worded:
Value: the value for this rating
#rating_data <- read.csv('OpioidMethRatingSubjectData.csv')
#clinical_data <- read.csv('DCR_ClinicalData.csv')
rating_data <- read.csv('DCR_RatingSubjectData-anonymized.csv')
clinical_data <- read.csv('DCR_Clinical-anonymized.csv')
subject_data <- merge(rating_data, clinical_data)library(ggplot2)
library(gridExtra)
library(grid)
library(ggpubr)## Warning: package 'ggpubr' was built under R version 3.5.3
## Loading required package: magrittr
#put levels in the order we like
subject_data$RelatedCategory <- factor(subject_data$RelatedCategory, levels = c('Neither', 'Meth', 'Opioid', 'Both'))
#subject_data$Category[subject_data$Category == 'meth_injection_instrument'] <- 'meth_instrument'
#cols <- c('control' = 'green', 'meth' = 'red', 'opioid' = 'blue', 'Control' = 'green', 'Meth' = 'red', 'Opioid' = 'blue')
cols <- c('control' = '#7CAE00', 'meth' = '#F8766D', 'opioid' = '#00BFC4', 'Control' = '#7CAE00', 'Meth' = '#F8766D', 'Opioid' = '#00BFC4')
for (subject in unique(subject_data$id)){
data_to_plot <- subject_data[subject_data$id == subject,]
p1 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'valence',]) +
geom_point(aes(x = Order, y = Value, color = ImageSet)) + ggtitle('Valence') +
xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) +
scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9)) +
theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
scale_color_manual(values = cols)
p2 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'arousal',]) +
geom_point(aes(x = Order, y = Value, color = ImageSet)) + ggtitle('Arousal') +
xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) +
scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9)) +
theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
scale_color_manual(values = cols)
p3 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'craving',]) +
geom_point(aes(x = Order, y = Value, color = ImageSet)) + ggtitle('Craving') +
xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) +
scale_y_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100)) +
theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
scale_color_manual(values = cols)
p4 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'typicality',]) +
geom_point(aes(x = Order, y = Value, color = ImageSet)) + ggtitle('Typicality') +
xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) +
scale_y_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100)) +
theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
scale_color_manual(values = cols)
p5 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'related',]) +
geom_point(aes(x = Order, y = RelatedCategory, color = ImageSet)) + ggtitle('Relatedness') +
xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) +
theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
scale_color_manual(values = cols)
filename <- paste(subject, 'responses.png', sep = '_')
png(filename, width = 1400, height = 400)
print(annotate_figure(ggarrange(p5, p1, p2, p3, p4, ncol = 5, nrow = 1, common.legend = TRUE, legend = 'right'), top = text_grob(subject, color = 'red')))
dev.off()
print(p1 + ggtitle(subject))
print(p2 + ggtitle(subject))
print(p3 + ggtitle(subject))
print(p4 + ggtitle(subject))
print(p5 + ggtitle(subject))
}library(reshape2)
one_summary <- function(dataset, prefix, grouping_field){
#dataset will be filtered down to either: all subjects, injection users, non-injection users
#prefix will be added to the output variable names, for example: allsubjects_arousal_mean or injectionusers_arousal_mean
#get raw mean for each image
mean_values <- aggregate(dataset[, c('Value')],
by = list(RatingType = dataset$RatingType, Group = dataset[, grouping_field]),
FUN = mean, na.rm = TRUE)
wide_mean_values <- dcast(mean_values, formula = Group ~ RatingType, value.var = 'x')
category_data <- dataset[dataset$RatingType == 'related',]
#so we can easily get a count of 'Neither', 'Both', 'Opioid', 'Meth'
category_data$count_col <- 1
#get counts of 'neither' ratings
counts_neither <- aggregate(category_data[category_data$RelatedCategory == 'Neither', c('count_col')],
by = list(Group = category_data[category_data$RelatedCategory == 'Neither', grouping_field]),
FUN = sum, na.rm = TRUE)
names(counts_neither) <- c('Group', 'CountNeither')
#'get count of 'meth' ratings
counts_meth <- aggregate(category_data[category_data$RelatedCategory == 'Meth', c('count_col')],
by = list(Group = category_data[category_data$RelatedCategory == 'Meth', grouping_field]),
FUN = sum, na.rm = TRUE)
names(counts_meth) <- c('Group', 'CountMeth')
all_counts <- merge(counts_neither, counts_meth, all = TRUE)
#get count of 'opioid' ratings
counts_opioid <- aggregate(category_data[category_data$RelatedCategory == 'Opioid', c('count_col')],
by = list(Group = category_data[category_data$RelatedCategory == 'Opioid', grouping_field]),
FUN = sum, na.rm = TRUE)
names(counts_opioid) <- c('Group', 'CountOpioid')
all_counts <- merge(all_counts, counts_opioid, all = TRUE)
#get count of 'both' ratings
counts_both <- aggregate(category_data[category_data$RelatedCategory == 'Both', c('count_col')],
by = list(Group = category_data[category_data$RelatedCategory == 'Both', grouping_field]),
FUN = sum, na.rm = TRUE)
names(counts_both) <- c('Group', 'CountBoth')
all_counts <- merge(all_counts, counts_both, all = TRUE)
#any 'NAs that appear here are ratings that never happened--e.g. control images that were never rated as related to meth
all_counts[is.na(all_counts)] <- 0
count_cols <- c('CountNeither', 'CountMeth', 'CountOpioid', 'CountBoth')
all_counts$FractionNeither <- all_counts$CountNeither / rowSums(all_counts[, count_cols])
all_counts$FractionMeth <- all_counts$CountMeth / rowSums(all_counts[, count_cols])
all_counts$FractionOpioid <- all_counts$CountOpioid / rowSums(all_counts[, count_cols])
all_counts$FractionBoth <- all_counts$CountBoth / rowSums(all_counts[, count_cols])
#this will be a number from -1 to 1: -1 for meth, 1 for opioid, and 0 for meth=opioid and many both/neither
all_counts$MethToOpioid <- (all_counts$CountOpioid - all_counts$CountMeth) / (all_counts$CountMeth + all_counts$CountOpioid + all_counts$CountBoth + all_counts$CountNeither)
#add prefix to all_counts column names
names(all_counts)[!(names(all_counts) %in% c('Group'))] <-
paste(prefix, names(all_counts)[!(names(all_counts) %in% c('Group'))], sep = '_')
#add _mean to the column names
names(wide_mean_values)[names(wide_mean_values) %in% c('arousal', 'craving', 'related', 'typicality', 'valence', 'valence_fromneutral')] <-
paste(prefix, names(wide_mean_values)[names(wide_mean_values) %in% c('arousal', 'craving', 'related', 'typicality', 'valence', 'valence_fromneutral')], 'mean', sep = '_')
#get raw SD for each image
sd_values <- aggregate(dataset[, c('Value')],
by = list(RatingType = dataset$RatingType, Group = dataset[, grouping_field]),
FUN = sd, na.rm = TRUE)
wide_sd_values <- dcast(sd_values, formula = Group ~ RatingType, value.var = 'x')
#add _sd to the column names
names(wide_sd_values)[names(wide_sd_values) %in% c('arousal', 'craving', 'related', 'typicality', 'valence', 'valence_fromneutral')] <-
paste(prefix, names(wide_sd_values)[names(wide_sd_values) %in% c('arousal', 'craving', 'related', 'typicality', 'valence', 'valence_fromneutral')], 'sd', sep = '_')
#merge them into one table
merged_summary_values <- merge(wide_mean_values, wide_sd_values)
merged_summary_values <- merge(merged_summary_values, all_counts)
return(merged_summary_values)
}Will report means/SDs for each image among all subjects, injection users, and non-injection users.
#also get abs(valence - 5) for each rating--distance from neutral
valence_rows <- subject_data[subject_data$RatingType == 'valence',]
valence_rows$RatingType <- 'valence_fromneutral'
valence_rows$Value <- abs(valence_rows$Value - 5)
subject_data <- rbind(subject_data, valence_rows)
all_data_summaries <- one_summary(subject_data, 'allsubjects', 'ImageSetFile')
injection_user_summaries <- one_summary(subject_data[subject_data$InjectionUser,], 'injectionusers', 'ImageSetFile')
noninjection_user_summaries <- one_summary(subject_data[!subject_data$InjectionUser,], 'noninjectionusers', 'ImageSetFile')
merged_summary_values <- merge(all_data_summaries, injection_user_summaries)
merged_summary_values <- merge(merged_summary_values, noninjection_user_summaries)
#now, add in HSV values computed by MATLAB
hsv_values <- read.csv('DCR_hsv_values.csv')
hsv_values$Group <- paste(hsv_values$ImageSet, hsv_values$File)
merged_summary_values <- merge(merged_summary_values, hsv_values)
#add in category column
category_data <- subject_data[!duplicated(subject_data$ImageSetFile), c('ImageSetFile', 'Category')]
names(category_data) <- c('Group', 'Category')
merged_summary_values <- merge(merged_summary_values, category_data)
merged_summary_values$p_cravingdelta_injectionVnon <- NA
###add in p-value for difference in craving between injection and non-injection users for each image
for (image in unique(subject_data$ImageSetFile)){
injection_craving_ratings <- subject_data[(subject_data$ImageSetFile == image) & subject_data$InjectionUser & subject_data$RatingType == 'craving', 'Value']
noninjection_craving_ratings <- subject_data[(subject_data$ImageSetFile == image) & (!subject_data$InjectionUser) & subject_data$RatingType == 'craving', 'Value']
this_test <- t.test(injection_craving_ratings, noninjection_craving_ratings)
merged_summary_values[merged_summary_values$Group == image, 'p_cravingdelta_injectionVnon'] <- this_test$p.value
}###add in category means/SDs for the same measures
all_data_summaries_categories <- one_summary(subject_data, 'allsubjects', 'Category')
injection_user_summaries_categories <- one_summary(subject_data[subject_data$InjectionUser,], 'injectionusers', 'Category')
noninjection_user_summaries_categories <- one_summary(subject_data[!subject_data$InjectionUser,], 'noninjectionusers', 'Category')
merged_summary_values_categories <- merge(all_data_summaries_categories, injection_user_summaries_categories)
merged_summary_values_categories <- merge(merged_summary_values_categories, noninjection_user_summaries_categories)
###add in mean/sd of hue/saturation/value mean/sd
this_mean_values <- aggregate(merged_summary_values[, c('hue_mean')],
by = list(Category = merged_summary_values$Category),
FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'hue_mean')
this_sd_values <- aggregate(merged_summary_values[, c('hue_mean')],
by = list(Category = merged_summary_values$Category),
FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'hue_sd')
merged_summary_values_categories <- merge(merged_summary_values_categories, this_mean_values)
merged_summary_values_categories <- merge(merged_summary_values_categories, this_sd_values)
##saturation
this_mean_values <- aggregate(merged_summary_values[, c('saturation_mean')],
by = list(Category = merged_summary_values$Category),
FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'saturation_mean')
this_sd_values <- aggregate(merged_summary_values[, c('saturation_mean')],
by = list(Category = merged_summary_values$Category),
FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'saturation_sd')
merged_summary_values_categories <- merge(merged_summary_values_categories, this_mean_values)
merged_summary_values_categories <- merge(merged_summary_values_categories, this_sd_values)
##value
this_mean_values <- aggregate(merged_summary_values[, c('value_mean')],
by = list(Category = merged_summary_values$Category),
FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'value_mean')
this_sd_values <- aggregate(merged_summary_values[, c('value_mean')],
by = list(Category = merged_summary_values$Category),
FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'value_sd')
merged_summary_values_categories <- merge(merged_summary_values_categories, this_mean_values)
merged_summary_values_categories <- merge(merged_summary_values_categories, this_sd_values)
#put ImageSet back in
library(stringr)
merged_summary_values_categories$ImageSet <- str_split_fixed(merged_summary_values_categories$Group, '_', 2)[,1]
###add in p-value for difference in craving between injection and non-injection users for each category
merged_summary_values_categories$p_cravingdelta_injectionVnon <- NA
for (category in unique(subject_data$Category)){
injection_craving_ratings <- subject_data[(subject_data$Category == category) & subject_data$InjectionUser & subject_data$RatingType == 'craving', 'Value']
noninjection_craving_ratings <- subject_data[(subject_data$Category == category) & (!subject_data$InjectionUser) & subject_data$RatingType == 'craving', 'Value']
this_test <- t.test(injection_craving_ratings, noninjection_craving_ratings)
merged_summary_values_categories[merged_summary_values_categories$Group == category, 'p_cravingdelta_injectionVnon'] <- this_test$p.value
}
library(dplyr)## Warning: package 'dplyr' was built under R version 3.5.3
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
##
## combine
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(data.table)##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## The following objects are masked from 'package:reshape2':
##
## dcast, melt
merged_summary_values <- rbindlist(list(merged_summary_values, merged_summary_values_categories), fill = TRUE, use.names = TRUE)
merged_summary_values <- data.frame(merged_summary_values)all_data_summaries_sets <- one_summary(subject_data, 'allsubjects', 'ImageSet')
injection_user_summaries_sets <- one_summary(subject_data[subject_data$InjectionUser,], 'injectionusers', 'ImageSet')
noninjection_user_summaries_sets <- one_summary(subject_data[!subject_data$InjectionUser,], 'noninjectionusers', 'ImageSet')
merged_summary_values_sets <- merge(all_data_summaries_sets, injection_user_summaries_sets)
merged_summary_values_sets <- merge(merged_summary_values_sets, noninjection_user_summaries_sets)
###HSV code
###add in mean/sd of hue/saturation/value mean/sd
this_mean_values <- aggregate(merged_summary_values[, c('hue_mean')],
by = list(ImageSet = merged_summary_values$ImageSet),
FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'hue_mean')
this_sd_values <- aggregate(merged_summary_values[, c('hue_mean')],
by = list(ImageSet = merged_summary_values$ImageSet),
FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'hue_sd')
merged_summary_values_sets <- merge(merged_summary_values_sets, this_mean_values)
merged_summary_values_sets <- merge(merged_summary_values_sets, this_sd_values)
##saturation
this_mean_values <- aggregate(merged_summary_values[, c('saturation_mean')],
by = list(ImageSet = merged_summary_values$ImageSet),
FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'saturation_mean')
this_sd_values <- aggregate(merged_summary_values[, c('saturation_mean')],
by = list(ImageSet = merged_summary_values$ImageSet),
FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'saturation_sd')
merged_summary_values_sets <- merge(merged_summary_values_sets, this_mean_values)
merged_summary_values_sets <- merge(merged_summary_values_sets, this_sd_values)
##value
this_mean_values <- aggregate(merged_summary_values[, c('value_mean')],
by = list(ImageSet = merged_summary_values$ImageSet),
FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'value_mean')
this_sd_values <- aggregate(merged_summary_values[, c('value_mean')],
by = list(ImageSet = merged_summary_values$ImageSet),
FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'value_sd')
merged_summary_values_sets <- merge(merged_summary_values_sets, this_mean_values)
merged_summary_values_sets <- merge(merged_summary_values_sets, this_sd_values)
#put ImageSet back in
library(stringr)
merged_summary_values_sets$ImageSet <- merged_summary_values_sets$Group
###add in p-value for difference in craving between injection and non-injection users for each category
merged_summary_values_sets$p_cravingdelta_injectionVnon <- NA
for (ImageSet in unique(subject_data$ImageSet)){
injection_craving_ratings <- subject_data[(subject_data$ImageSet == ImageSet) & subject_data$InjectionUser & subject_data$RatingType == 'craving', 'Value']
noninjection_craving_ratings <- subject_data[(subject_data$ImageSet == ImageSet) & (!subject_data$InjectionUser) & subject_data$RatingType == 'craving', 'Value']
this_test <- t.test(injection_craving_ratings, noninjection_craving_ratings)
merged_summary_values_sets[merged_summary_values_sets$Group == ImageSet, 'p_cravingdelta_injectionVnon'] <- this_test$p.value
}
library(dplyr)
library(data.table)
merged_summary_values <- rbindlist(list(merged_summary_values, merged_summary_values_sets), fill = TRUE, use.names = TRUE)
merged_summary_values <- data.frame(merged_summary_values)
###HSV code#variables to write out, in a reasonable order
vars_to_save <- c("Group", "ImageSet", "File", "Category",
"allsubjects_valence_mean", "allsubjects_valence_sd", "allsubjects_arousal_mean", "allsubjects_arousal_sd",
"allsubjects_craving_mean", "allsubjects_craving_sd", "allsubjects_typicality_mean", "allsubjects_typicality_sd",
"allsubjects_FractionNeither", "allsubjects_FractionMeth", "allsubjects_FractionOpioid", "allsubjects_FractionBoth", "allsubjects_MethToOpioid",
"injectionusers_valence_mean", "injectionusers_valence_sd", "injectionusers_arousal_mean", "injectionusers_arousal_sd",
"injectionusers_craving_mean", "injectionusers_craving_sd", "injectionusers_typicality_mean", "injectionusers_typicality_sd",
"injectionusers_FractionNeither", "injectionusers_FractionMeth", "injectionusers_FractionOpioid", "injectionusers_FractionBoth", "injectionusers_MethToOpioid",
"noninjectionusers_valence_mean", "noninjectionusers_valence_sd", "noninjectionusers_arousal_mean", "noninjectionusers_arousal_sd",
"noninjectionusers_craving_mean", "noninjectionusers_craving_sd", "noninjectionusers_typicality_mean", "noninjectionusers_typicality_sd",
"noninjectionusers_FractionNeither", "noninjectionusers_FractionMeth", "noninjectionusers_FractionOpioid", "noninjectionusers_FractionBoth",
"noninjectionusers_MethToOpioid",
"hue_mean", "hue_sd", "saturation_mean", "saturation_sd", "value_mean", "value_sd", "p_cravingdelta_injectionVnon")
write.csv(format(merged_summary_values[, vars_to_save], digits = 3), 'DCR_summaries_7-3-2019.csv', row.names = FALSE)
#only want image summaries, not category means, for example
plot_vals <- merged_summary_values[!is.na(merged_summary_values$File),]
#see how distance from neutral compares with mean valence
ggplot(data = plot_vals, aes(x = allsubjects_valence_mean, y = allsubjects_valence_fromneutral_mean, color = ImageSet)) + geom_point() +
xlab('valence') + ylab('valence distance to mean')+ geom_smooth(method = 'lm') +
scale_color_manual(values = cols)print("#####Test for control images")## [1] "#####Test for control images"
cor.test(plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == 'control'],plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == 'control'] )##
## Pearson's product-moment correlation
##
## data: plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "control"] and plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "control"] and "control"]
## t = -0.0043273, df = 118, p-value = 0.9966
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1796269 0.1788557
## sample estimates:
## cor
## -0.0003983601
print("#####Test for meth images")## [1] "#####Test for meth images"
cor.test(plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == 'meth'],plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == 'meth'] )##
## Pearson's product-moment correlation
##
## data: plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "meth"] and plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "meth"] and "meth"]
## t = -6.3426, df = 118, p-value = 4.322e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6268160 -0.3572733
## sample estimates:
## cor
## -0.504225
print("#####Test for opioid images")## [1] "#####Test for opioid images"
cor.test(plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == 'opioid'],plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == 'opioid'] )##
## Pearson's product-moment correlation
##
## data: plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "opioid"] and plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "opioid"] and "opioid"]
## t = -11.6, df = 118, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8039688 -0.6335581
## sample estimates:
## cor
## -0.7299109
scatter_data <- subject_data[subject_data$RatingType == c('valence'), c('id', 'ImageSet', 'ImageSetFile', 'RatingType', 'Value')]
names(scatter_data)[names(scatter_data) == 'Value'] <- 'Valence'
scatter_data$ValenceFromMean <- abs(scatter_data$Valence - 5)
#c <- cor(craving_data$Craving, typicality_data$Typicality)
ggplot(scatter_data, aes(x = Valence, y = ValenceFromMean, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm') +
scale_color_manual(values = cols)ggplot(scatter_data, aes(x = Valence, y = ValenceFromMean, color = ImageSet)) + geom_hex() + geom_smooth(method = 'lm') +
scale_color_manual(values = cols)## Warning: package 'hexbin' was built under R version 3.5.3
scatter_data <- subject_data[subject_data$ImageSet %in% c('meth', 'opioid'), c('id', 'ImageSet', 'ImageSetFile', 'RatingType', 'Value')]
craving_data <- scatter_data[scatter_data$RatingType == 'craving',]
typicality_data <- scatter_data[scatter_data$RatingType == 'typicality',]
craving_data$Craving <- craving_data$Value
typicality_data$Typicality <- typicality_data$Value
correlation <- cor(craving_data$Craving, typicality_data$Typicality)
scatter_data <- merge(craving_data[, c('id', 'ImageSet', 'ImageSetFile', 'Craving')], typicality_data[, c('id', 'ImageSet', 'ImageSetFile', 'Typicality')])
ggplot(scatter_data, aes(x = Craving, y = Typicality, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm') +
ggtitle(paste0('r = ', round(correlation, 2), '\nBut not appropriate, since multiple measures per subject')) +
scale_color_manual(values = cols)mean_scatter_data <- merged_summary_values[merged_summary_values$ImageSet %in% c('meth', 'opioid'),]
correlation <- cor(mean_scatter_data$allsubjects_craving_mean, mean_scatter_data$allsubjects_typicality_mean)
ggplot(mean_scatter_data, aes(x = allsubjects_craving_mean, y = allsubjects_typicality_mean, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm') +
scale_y_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100)) +
scale_x_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100)) +
xlab('Craving Mean') + ylab ('Typicality Mean') + ggtitle(paste0('r = ', round(correlation, 2))) +
scale_color_manual(values = cols)make_plotset <- function(data_to_plot, main_title, filename, group){
p1 <- ggplot(data_to_plot, aes_string(x = 'allsubjects_valence_mean', fill = 'ImageSet')) +
geom_histogram(binwidth = 0.5, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
scale_color_manual(values = cols)
p2 <- ggplot(data_to_plot, aes_string(x = 'allsubjects_arousal_mean', fill = 'ImageSet')) +
geom_histogram(binwidth = 0.5, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
scale_color_manual(values = cols)
p3 <- ggplot(data_to_plot, aes_string(x = 'allsubjects_craving_mean', fill = 'ImageSet')) +
geom_histogram(binwidth = 10, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
scale_color_manual(values = cols)
p4 <- ggplot(data_to_plot, aes_string(x = 'allsubjects_typicality_mean', fill = 'ImageSet')) +
geom_histogram(binwidth = 10, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
scale_color_manual(values = cols)
p5 <- ggplot(data_to_plot, aes_string(x = 'hue_mean', fill = 'ImageSet')) +
geom_histogram(binwidth = 0.1, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
scale_color_manual(values = cols)
p6 <- ggplot(data_to_plot, aes_string(x = 'saturation_mean', fill = 'ImageSet')) +
geom_histogram(binwidth = 0.1, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
scale_color_manual(values = cols)
p7 <- ggplot(data_to_plot, aes_string(x = 'value_mean', fill = 'ImageSet')) +
geom_histogram(binwidth = 0.1, alpha = 0.5, position = 'identity') +
scale_color_manual(values = cols)
#png(filename, width = 1400, height = 400)
#grid.arrange(p1, p2, p3, p4, p5, p6, p7, ncol = 7, top = textGrob(main_title, gp=gpar(fontsize = 50)))
#dev.off()
print(p1)
print(p2)
print(p3)
print(p4)
print(p5)
print(p6)
print(p7)
}
library(sm)## Warning: package 'sm' was built under R version 3.5.3
## Package 'sm', version 2.2-5.6: type help(sm) for summary information
one_density_plot <- function(var_to_plot, label){
sm.density.compare(merged_summary_values[, var_to_plot], as.factor(merged_summary_values$ImageSet), lwd = 5)
title(main=label)
colfill<-c(2:(2+length(levels(as.factor(merged_summary_values$ImageSet)))))
legend('topright', levels(as.factor(merged_summary_values$ImageSet)), fill=colfill)
}
source('R_rainclouds.R')
library(cowplot)## Warning: package 'cowplot' was built under R version 3.5.3
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
##
## get_legend
## The following object is masked from 'package:ggplot2':
##
## ggsave
one_raincloud_plot <- function(var_to_plot, label){
p6 <- ggplot(subject_data[subject_data$RatingType == 'arousal',],aes(x=ImageSet,y=Value, fill = ImageSet, colour = ImageSet))+
geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim = FALSE)+
#geom_point(position = position_jitter(width = .15), size = .25)+
geom_point(position = position_jitter(width = 1), size = .25)+
geom_boxplot(aes(x = as.numeric(ImageSet)+0.25, y = Value),outlier.shape = NA,
alpha = 0.3, width = .1, colour = "BLACK") +
ylab('Arousal')+xlab('ImageSet')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
colour = FALSE) +
ggtitle("Figure R6: Change in Colour Palette") +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols)
#scale_colour_brewer(palette = "Dark2")+
#scale_fill_brewer(palette = "Dark2")+
ggsave(paste0(label, '.png'), width = w, height = h)
p6
}
subject_data$ImageSetCap <- as.character(subject_data$ImageSet)
subject_data$ImageSetCap[subject_data$ImageSetCap == 'meth'] <- 'Meth'
subject_data$ImageSetCap[subject_data$ImageSetCap == 'opioid'] <- 'Opioid'
subject_data$ImageSetCap[subject_data$ImageSetCap == 'control'] <- 'Control'
subject_data$ImageSetCap <- factor(subject_data$ImageSetCap, levels = c('Opioid', 'Meth', 'Control'))
p1 <- ggplot(subject_data[subject_data$RatingType == 'craving',],aes(x=ImageSetCap,y=Value, fill = ImageSetCap, colour = ImageSetCap))+
geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim = FALSE)+
geom_point(position = position_jitter(width = .15), size = .25)+
geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = Value),outlier.shape = NA,
alpha = 0.3, width = .1, colour = "BLACK") +
ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
colour = FALSE) +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols) +
ggtitle("Craving") + scale_y_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100))
p2 <- ggplot(subject_data[subject_data$RatingType == 'arousal',],aes(x=ImageSetCap,y=Value, fill = ImageSetCap, colour = ImageSetCap))+
geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim = FALSE)+
geom_point(position = position_jitter(width = .15), size = .25)+
geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = Value),outlier.shape = NA,
alpha = 0.3, width = .1, colour = "BLACK") +
ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
colour = FALSE) +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols) +
ggtitle("Arousal") + scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))
p3 <- ggplot(subject_data[subject_data$RatingType == 'valence',],aes(x=ImageSetCap,y=Value, fill = ImageSetCap, colour = ImageSetCap))+
geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim = FALSE)+
geom_point(position = position_jitter(width = .15), size = .25)+
geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = Value),outlier.shape = NA,
alpha = 0.3, width = .1, colour = "BLACK") +
ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
colour = FALSE) +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols) +
ggtitle("Valence") + scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))
png('rainclouds.png', width = 1400, height = 400)
print(annotate_figure(ggarrange(p1, p2, p3, ncol = 3, nrow = 1, common.legend = TRUE, legend = 'right'), top = text_grob('Rainclouds', color = 'red')))## Warning: Removed 168 rows containing missing values (geom_flat_violin).
## Warning: Removed 1858 rows containing missing values (geom_point).
## Warning: Removed 380 rows containing missing values (geom_flat_violin).
## Warning: Removed 2096 rows containing missing values (geom_point).
## Warning: Removed 252 rows containing missing values (geom_flat_violin).
## Warning: Removed 464 rows containing missing values (geom_point).
## Warning: Removed 168 rows containing missing values (geom_flat_violin).
## Warning: Removed 1858 rows containing missing values (geom_point).
## Warning: Removed 380 rows containing missing values (geom_flat_violin).
## Warning: Removed 2096 rows containing missing values (geom_point).
## Warning: Removed 252 rows containing missing values (geom_flat_violin).
## Warning: Removed 464 rows containing missing values (geom_point).
dev.off()## png
## 2
p4 <- ggplot(subject_data,aes(x=ImageSetCap,y=millisecondsOnPage / 1000, fill = ImageSetCap, colour = ImageSetCap))+
geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim = FALSE)+
geom_point(position = position_jitter(width = .15), size = .25)+
geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = millisecondsOnPage / 1000),outlier.shape = NA,
alpha = 0.3, width = .1, colour = "BLACK") +
ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
colour = FALSE) +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols) +
ggtitle("SecondsOnPage") #+ scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))
p5 <- ggplot(subject_data,aes(x=ImageSetCap,y=responseDelay / 1000, fill = ImageSetCap, colour = ImageSetCap))+
geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim = FALSE)+
geom_point(position = position_jitter(width = .15), size = .25)+
geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = responseDelay / 1000),outlier.shape = NA,
alpha = 0.3, width = .1, colour = "BLACK") +
ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
colour = FALSE) +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols) +
ggtitle("ResponseDelay") #+ scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))
png('RTrainclouds.png', width = 1400, height = 400)
print(annotate_figure(ggarrange(p4, p5, ncol = 1, nrow = 2, common.legend = TRUE, legend = 'right'), top = text_grob('RTRainclouds', color = 'red')))
dev.off()## png
## 2
p1 <- ggplot(subject_data[subject_data$RatingType == 'craving',],
aes(x = Value, y = millisecondsOnPage / 1000,color = ImageSet)) + geom_point() +
xlab('Craving') + ylab('ScondsOnPage') + #scale_x_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) +
#scale_y_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) +
theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
geom_smooth(method = 'lm') +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols)
p2 <- ggplot(subject_data[subject_data$RatingType == 'craving',],
aes(x = Value, y = responseDelay / 1000,color = ImageSet)) + geom_point() +
xlab('Craving') + ylab('ResponseDelay') + #scale_x_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) +
#scale_y_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) +
theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
geom_smooth(method = 'lm') +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols)
png('RTCravingScatter.png', width = 400, height = 800)
print(annotate_figure(ggarrange(p1, p2, ncol = 1, nrow = 2, common.legend = TRUE, legend = 'right'), top = text_grob('RTScatters', color = 'red')))
dev.off()## png
## 2
craving_data <- subject_data[subject_data$RatingType == 'craving',]
#set threshold of 20 seconds, remove 0.6% of data
sum(craving_data$millisecondsOnPage > 10000) / nrow(craving_data)## [1] 0.02434343
one_density_plot('allsubjects_valence_mean', 'Valence')one_density_plot('allsubjects_arousal_mean', 'Arousal')one_density_plot('allsubjects_craving_mean', 'Craving')one_density_plot('allsubjects_typicality_mean', 'Typicality')one_density_plot('hue_mean', 'Hue')one_density_plot('saturation_mean', 'Saturation')one_density_plot('value_mean', 'Value')library(corrplot)## corrplot 0.84 loaded
library(PerformanceAnalytics)## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Attaching package: 'xts'
## The following objects are masked from 'package:data.table':
##
## first, last
## The following objects are masked from 'package:dplyr':
##
## first, last
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
cols_to_use <- c('allsubjects_valence_mean', 'allsubjects_valence_fromneutral_mean', 'allsubjects_arousal_mean',
'allsubjects_craving_mean')
one_corrplot <- function(dset, title){
col2 <- colorRampPalette(rev(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
"#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
"#4393C3", "#2166AC", "#053061")))
corvals <- cor(dset)
rownames(corvals) <- colnames(corvals) <- c('valence', 'abs_valence', 'arousal', 'craving')
corrplot.mixed(corvals, upper = 'ellipse', lower = 'number', title = title, mar=c(0,0,2,0),
upper.col = col2(50), lower.col = col2(50))
chart.Correlation(dset)
#returns z-transformed correlations, for easier comparisons between matrices
return(atanh(corvals))
}
p <- ggplot(merged_summary_values[merged_summary_values$ImageSet %in% c('meth', 'opioid'),], aes(x = allsubjects_craving_mean, y = allsubjects_typicality_mean,
color = ImageSet)) + geom_point() +
xlab('Craving') + ylab('Typicality') + scale_x_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) +
scale_y_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) +
theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
geom_smooth(method = 'lm') +
scale_color_manual(values = cols)
#wide_subject_data <- dcast(subject_data, formula = id ~ RatingType + ImageSetFile, value.var = 'Value')
p <- ggplot(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),],
aes(x = allsubjects_craving_mean, y = allsubjects_typicality_mean,color = ImageSet)) + geom_point() +
xlab('Craving') + ylab('Typicality') + scale_x_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) +
scale_y_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) +
theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
geom_smooth(method = 'lm') +
scale_color_manual(values = cols)
drug_zs <- one_corrplot(merged_summary_values[merged_summary_values$ImageSet %in% c('meth', 'opioid'), cols_to_use], 'Meth/Opioid')meth_zs <- one_corrplot(merged_summary_values[merged_summary_values$ImageSet %in% c('meth'), cols_to_use], 'Meth')opioid_zs <- one_corrplot(merged_summary_values[merged_summary_values$ImageSet %in% c('opioid'), cols_to_use], 'Opioid')control_zs <- one_corrplot(merged_summary_values[merged_summary_values$ImageSet %in% c('control'), cols_to_use], 'control')#to compare correlation z-scores--here we use N1=N2=120, the number of images
#Zobserved = (z1 - z2) / (square root of [ (1 / N1 - 3) + (1 / N2 - 3) ])
meth_opioid_deltazs <- (meth_zs - opioid_zs) / sqrt(1/(120-3) + 1/(120-3))
meth_opioid_deltaps <- 2*pnorm(-abs(meth_opioid_deltazs))
col2 <- colorRampPalette(rev(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
"#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
"#4393C3", "#2166AC", "#053061")))
corrplot.mixed(meth_opioid_deltazs, upper = 'color', lower = 'number', title = 'Delta Z Score (Meth - Opioid)', mar=c(0,0,2,0), is.corr = FALSE,
p.mat = meth_opioid_deltaps, sig.level = 1, lower.col = col2(50), upper.col = col2(50))corrplot.mixed(meth_opioid_deltaps, upper = 'color', lower = 'number', title = 'Delta r p-values (Meth - Opioid)', mar=c(0,0,2,0), is.corr = FALSE,
lower.col = col2(50), upper.col = col2(50))# p.mat = meth_opioid_deltaps, sig.level = 0.01, lower.col = col2(50), upper.col = col2(50))
print('Matrix of p-values for Correlation Differences')## [1] "Matrix of p-values for Correlation Differences"
print(meth_opioid_deltaps)## valence abs_valence arousal craving
## valence NaN 0.003167665 0.005306939 0.0232624
## abs_valence 0.003167665 NaN 0.130490960 0.7358989
## arousal 0.005306939 0.130490960 NaN 0.2669365
## craving 0.023262397 0.735898874 0.266936547 NaN
print('Matrix of fdr corrected p-values')## [1] "Matrix of fdr corrected p-values"
p.adjust(meth_opioid_deltaps[upper.tri(meth_opioid_deltaps)], method = 'fdr')## [1] 0.01592082 0.01592082 0.19573644 0.04652479 0.73589887 0.32032386
meth_craving ~ Age + sex + age_of_meth_onset + duration_of_meth_use + cost_of_meth + days_of_meth + duration_of_current_abstinance
#first, corrplot of these variables
clinical_columns <- c('Age', 'SexMale', 'InjectionUser', 'MethOnsetAge', 'MethMonthsUsed', 'MethDollars',
'MethDaysLastMonth','MethDaysAbstinant', 'OpioidOnsetAge', 'OpioidMonthsUsed',
'OpioidDaysLastMonth', 'OpioidDaysAbstinant', 'HeroinOnsetAge', 'HeroinMonthsUsed',
'HeroinDaysLastMonth', 'HeroinDaysAbstinant', 'OpioidOrHeroinOnsetAge',
'OpioidOrHeroinMonthsUsed', 'OpioidOrHeroinDollars', 'OpioidOrHeroinDaysLastMonth','OpioidOrHeroinDaysAbstinant')
col2 <- colorRampPalette(rev(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
"#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
"#4393C3", "#2166AC", "#053061")))
clinical_data$SexMale <- 0
clinical_data$SexMale[clinical_data$Sex == 'Male'] <- 1
corvals <- cor(clinical_data[, clinical_columns], use = 'pairwise.complete')
png('ClinicalCorrelations.png', width = 1000, height = 1000)
corrplot.mixed(corvals, upper = 'ellipse', lower = 'number', title = 'Clinical Variable Correlations', mar=c(0,0,2,0),
upper.col = col2(50), lower.col = col2(50), tl.pos = 'lt')
dev.off()## png
## 2
png('ClinicalScatterplotMatrix.png', width = 1000, height = 1000)
chart.Correlation(clinical_data[, clinical_columns])
dev.off()## png
## 2
library(tableone)
CreateTableOne(vars = clinical_columns, data = clinical_data)##
## Overall
## n 28
## Age (mean (sd)) 37.71 (8.11)
## SexMale (mean (sd)) 0.57 (0.50)
## InjectionUser = TRUE (%) 21 (75.0)
## MethOnsetAge (mean (sd)) 19.42 (5.40)
## MethMonthsUsed (mean (sd)) 117.96 (96.28)
## MethDollars (mean (sd)) 2928.93 (2540.84)
## MethDaysLastMonth (mean (sd)) 21.00 (12.29)
## MethDaysAbstinant (mean (sd)) 729.32 (858.29)
## OpioidOnsetAge (mean (sd)) 20.82 (6.39)
## OpioidMonthsUsed (mean (sd)) 115.32 (94.96)
## OpioidDaysLastMonth (mean (sd)) 16.71 (12.91)
## OpioidDaysAbstinant (mean (sd)) 751.04 (664.18)
## HeroinOnsetAge (mean (sd)) 26.79 (7.73)
## HeroinMonthsUsed (mean (sd)) 31.43 (71.31)
## HeroinDaysLastMonth (mean (sd)) 10.61 (13.42)
## HeroinDaysAbstinant (mean (sd)) 1492.20 (2022.00)
## OpioidOrHeroinOnsetAge (mean (sd)) 20.54 (6.29)
## OpioidOrHeroinMonthsUsed (mean (sd)) 117.89 (92.67)
## OpioidOrHeroinDollars (mean (sd)) 1773.21 (1711.63)
## OpioidOrHeroinDaysLastMonth (mean (sd)) 22.36 (11.18)
## OpioidOrHeroinDaysAbstinant (mean (sd)) 914.43 (1332.91)
library(nlme)##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
library(sjmisc)## Warning: package 'sjmisc' was built under R version 3.5.3
library(sjPlot)## Warning: package 'sjPlot' was built under R version 3.5.3
##
## Attaching package: 'sjPlot'
## The following objects are masked from 'package:cowplot':
##
## plot_grid, save_plot
library(lme4)## Loading required package: Matrix
##
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
##
## lmList
library(lmerTest)## Warning: package 'lmerTest' was built under R version 3.5.3
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
this_lme <- lmer(Value ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit +(1 | id),
data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'meth',])## Warning: Some predictor variables are on very different scales: consider rescaling
## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'meth craving', show.values = TRUE, type = 'std', order.terms = c(1,2,8,9,10,3,4,6,5,7)))print(summary(this_lme))## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars +
## MethDaysLastMonth + MethDaysAbstinant + Time * Visit + (1 | id)
## Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==
## "meth", ]
##
## REML criterion at convergence: 26081.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.9873 -0.1122 0.1452 0.5159 2.9383
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 132.9 11.53
## Residual 204.4 14.30
## Number of obs: 3180, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.467e+01 1.915e+01 1.907e+01 4.421 0.000291 ***
## Age 4.888e-01 3.550e-01 1.901e+01 1.377 0.184622
## SexMale -8.208e+00 5.911e+00 1.901e+01 -1.389 0.181023
## MethOnsetAge 7.246e-02 5.766e-01 1.902e+01 0.126 0.901319
## MethMonthsUsed 9.298e-03 3.612e-02 1.900e+01 0.257 0.799619
## MethDollars 1.451e-03 1.191e-03 1.904e+01 1.218 0.238067
## MethDaysLastMonth -5.150e-01 2.326e-01 1.902e+01 -2.213 0.039282 *
## MethDaysAbstinant -4.520e-03 8.542e-03 1.900e+01 -0.529 0.602841
## Time -2.011e-02 6.679e-03 3.150e+03 -3.011 0.002625 **
## VisitV2 3.247e-01 1.029e+00 3.151e+03 0.316 0.752341
## Time:VisitV2 -4.525e-02 9.591e-03 3.150e+03 -4.717 2.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Age SexMal MthOnA MthMnU MthDll MthDLM MthDyA Time VistV2
## Age -0.526
## SexMale -0.061 -0.002
## MethOnsetAg -0.570 -0.154 -0.386
## MthMnthsUsd -0.222 -0.384 -0.330 0.580
## MethDollars -0.337 0.114 -0.009 0.166 -0.053
## MthDysLstMn -0.155 -0.110 0.346 -0.016 -0.074 -0.433
## MthDysAbstn -0.507 0.072 0.334 0.051 0.101 0.386 0.058
## Time -0.032 0.000 -0.002 0.002 0.002 -0.002 -0.001 -0.004
## VisitV2 -0.022 -0.004 0.002 -0.002 0.003 -0.003 0.001 -0.003 0.602
## Time:VistV2 0.021 0.002 0.000 0.000 -0.002 0.000 0.002 0.002 -0.695 -0.868
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Age 387.3 387.3 1 19.01 1.8951 0.18462
## Sex 394.0 394.0 1 19.01 1.9281 0.18102
## MethOnsetAge 3.2 3.2 1 19.02 0.0158 0.90132
## MethMonthsUsed 13.5 13.5 1 19.00 0.0663 0.79962
## MethDollars 303.2 303.2 1 19.04 1.4838 0.23807
## MethDaysLastMonth 1001.3 1001.3 1 19.02 4.8996 0.03928 *
## MethDaysAbstinant 57.2 57.2 1 19.00 0.2800 0.60284
## Time 16175.7 16175.7 1 3150.40 79.1509 < 2.2e-16 ***
## Visit 20.4 20.4 1 3150.88 0.0996 0.75234
## Time:Visit 4547.9 4547.9 1 3150.24 22.2536 2.493e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
scaled_data <- subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'meth',]
scaled_data[, c('Age', 'MethOnsetAge', 'MethMonthsUsed', 'MethDollars', 'MethDaysLastMonth', 'MethDaysAbstinant')] <-
scale(scaled_data[, c('Age', 'MethOnsetAge', 'MethMonthsUsed', 'MethDollars', 'MethDaysLastMonth', 'MethDaysAbstinant')])
this_lme <- lmer(Value ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit + (1 | id),
data = scaled_data)
print(plot_model(this_lme, title = 'meth craving scaled predictors', show.values = TRUE, type = 'std', order.terms = c(1,2,8,9,10,3,4,6,5,7)))print(summary(this_lme))## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars +
## MethDaysLastMonth + MethDaysAbstinant + Time * Visit + (1 | id)
## Data: scaled_data
##
## REML criterion at convergence: 26030.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.9873 -0.1122 0.1452 0.5159 2.9383
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 132.9 11.53
## Residual 204.4 14.30
## Number of obs: 3180, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.582e+01 4.158e+00 2.003e+01 23.041 6.85e-16 ***
## Age 3.926e+00 2.852e+00 1.901e+01 1.377 0.18462
## SexMale -8.208e+00 5.911e+00 1.901e+01 -1.389 0.18102
## MethOnsetAge 3.795e-01 3.020e+00 1.902e+01 0.126 0.90132
## MethMonthsUsed 8.708e-01 3.383e+00 1.900e+01 0.257 0.79962
## MethDollars 3.651e+00 2.998e+00 1.904e+01 1.218 0.23807
## MethDaysLastMonth -6.240e+00 2.819e+00 1.902e+01 -2.213 0.03928 *
## MethDaysAbstinant -3.834e+00 7.246e+00 1.900e+01 -0.529 0.60284
## Time -2.011e-02 6.679e-03 3.150e+03 -3.011 0.00262 **
## VisitV2 3.247e-01 1.029e+00 3.151e+03 0.316 0.75234
## Time:VisitV2 -4.525e-02 9.591e-03 3.150e+03 -4.717 2.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Age SexMal MthOnA MthMnU MthDll MthDLM MthDyA Time VistV2
## Age 0.076
## SexMale -0.761 -0.002
## MethOnsetAg 0.371 -0.154 -0.386
## MthMnthsUsd 0.336 -0.384 -0.330 0.580
## MethDollars 0.136 0.114 -0.009 0.166 -0.053
## MthDysLstMn -0.298 -0.110 0.346 -0.016 -0.074 -0.433
## MthDysAbstn 0.047 0.072 0.334 0.051 0.101 0.386 0.058
## Time -0.148 0.000 -0.002 0.002 0.002 -0.002 -0.001 -0.004
## VisitV2 -0.122 -0.004 0.002 -0.002 0.003 -0.003 0.001 -0.003 0.602
## Time:VistV2 0.104 0.002 0.000 0.000 -0.002 0.000 0.002 0.002 -0.695 -0.868
print(anova(this_lme))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Age 387.3 387.3 1 19.01 1.8951 0.18462
## Sex 394.0 394.0 1 19.01 1.9281 0.18102
## MethOnsetAge 3.2 3.2 1 19.02 0.0158 0.90132
## MethMonthsUsed 13.5 13.5 1 19.00 0.0663 0.79962
## MethDollars 303.2 303.2 1 19.04 1.4838 0.23807
## MethDaysLastMonth 1001.3 1001.3 1 19.02 4.8996 0.03928 *
## MethDaysAbstinant 57.2 57.2 1 19.00 0.2800 0.60284
## Time 16175.7 16175.7 1 3150.40 79.1509 < 2.2e-16 ***
## Visit 20.4 20.4 1 3150.88 0.0996 0.75234
## Time:Visit 4547.9 4547.9 1 3150.24 22.2536 2.493e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
this_lme <- lmer(Value ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed + OpioidOrHeroinDollars +
OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant + Time * Visit + (1 | id),
data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'opioid',])## Warning: Some predictor variables are on very different scales: consider rescaling
## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'opioid craving', show.values = TRUE, type = 'std', order.terms = c(1,2,8, 9, 10, 3,4,6,5,7)))print(summary(this_lme))## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed +
## OpioidOrHeroinDollars + OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant +
## Time * Visit + (1 | id)
## Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==
## "opioid", ]
##
## REML criterion at convergence: 25942.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3037 -0.1131 0.1809 0.4864 3.1743
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 139.1 11.80
## Residual 195.4 13.98
## Number of obs: 3180, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.189e+01 1.626e+01 1.905e+01 5.036 7.28e-05 ***
## Age 2.123e-01 4.675e-01 1.900e+01 0.454 0.654880
## SexMale -5.248e+00 5.059e+00 1.903e+01 -1.037 0.312565
## OpioidOrHeroinOnsetAge 1.731e-01 5.309e-01 1.904e+01 0.326 0.747881
## OpioidOrHeroinMonthsUsed 5.140e-02 3.799e-02 1.902e+01 1.353 0.191893
## OpioidOrHeroinDollars 1.204e-03 1.784e-03 1.904e+01 0.675 0.507923
## OpioidOrHeroinDaysLastMonth -2.161e-01 2.839e-01 1.900e+01 -0.761 0.455866
## OpioidOrHeroinDaysAbstinant 3.980e-04 2.047e-03 1.901e+01 0.194 0.847885
## Time -3.784e-02 6.548e-03 3.150e+03 -5.779 8.27e-09 ***
## VisitV2 -3.337e+00 9.908e-01 3.151e+03 -3.368 0.000767 ***
## Time:VisitV2 -1.299e-02 9.328e-03 3.150e+03 -1.392 0.164007
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Age SexMal OpOHOA OpOHMU OpdOHD OOHDLM OpOHDA Time VistV2
## Age -0.756
## SexMale -0.067 0.019
## OpdOrHrnOnA 0.085 -0.616 -0.142
## OpdOrHrnMnU 0.128 -0.396 -0.332 0.405
## OpdOrHrnDll 0.057 0.117 0.209 -0.417 -0.245
## OpdOrHrnDLM -0.633 0.487 0.057 -0.264 -0.393 -0.157
## OpdOrHrnDyA -0.177 -0.044 -0.157 -0.007 0.225 0.020 0.235
## Time -0.034 -0.001 0.001 -0.002 -0.001 0.001 -0.002 0.000
## VisitV2 -0.026 -0.001 0.002 -0.003 -0.003 0.004 -0.002 -0.001 0.603
## Time:VistV2 0.023 0.003 0.001 -0.001 0.001 -0.001 0.002 -0.001 -0.702 -0.863
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Age 40.3 40.3 1 19.00 0.2062 0.6548801
## Sex 210.3 210.3 1 19.03 1.0761 0.3125645
## OpioidOrHeroinOnsetAge 20.8 20.8 1 19.04 0.1064 0.7478808
## OpioidOrHeroinMonthsUsed 357.7 357.7 1 19.02 1.8309 0.1918931
## OpioidOrHeroinDollars 89.0 89.0 1 19.04 0.4553 0.5079235
## OpioidOrHeroinDaysLastMonth 113.2 113.2 1 19.00 0.5795 0.4558659
## OpioidOrHeroinDaysAbstinant 7.4 7.4 1 19.01 0.0378 0.8478851
## Time 17646.5 17646.5 1 3150.25 90.3210 < 2.2e-16 ***
## Visit 2215.8 2215.8 1 3150.89 11.3411 0.0007672 ***
## Time:Visit 378.6 378.6 1 3150.27 1.9378 0.1640067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We will look drug and neutral images separately. For drug images, the model will be:
Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)
And for neutral images it will be:
Value ~ Time * Visit + (Age * Sex) + (1 | id)
for valence, arousal, and craving separately
summarize_one_lme_drug <- function(dset, RatingType){
this_lme <- lmer(Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id),
data = dset[dset$RatingType == RatingType,])
print(plot_model(this_lme, title = RatingType, show.values = TRUE, type = 'std'))
print(summary(this_lme))
#for using nlme instead of lme4--but can't get standardized betas in the forest plot
#that_lme <- lme(fixed = Value ~ Time * Visit + Age + Sex + ImageSet, random = ~ 1 | id,
# data = dset[dset$RatingType == RatingType,])
#print(summary(that_lme))
print(anova(this_lme))
}
summarize_one_lme_neutral <- function(dset, RatingType){
this_lme <- lmer(Value ~ Time * Visit + (Age * Sex) + (1 | id),
data = dset[dset$RatingType == RatingType,])
print(plot_model(this_lme, title = RatingType, show.values = TRUE, type = 'std'))
print(summary(this_lme))
print(anova(this_lme))
}summarize_one_lme_drug(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'valence')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 18920.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8998 -0.3812 0.0481 0.4229 6.6748
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 2.041 1.4287
## Residual 0.995 0.9975
## Number of obs: 6597, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.863e+00 1.687e+00 2.412e+01 2.290 0.0311 *
## Time 2.065e-03 3.230e-04 6.562e+03 6.393 1.74e-10 ***
## VisitV2 3.917e-01 4.943e-02 6.562e+03 7.925 2.66e-15 ***
## Age -3.369e-03 4.350e-02 2.410e+01 -0.077 0.9389
## SexMale -5.422e+00 2.682e+00 2.410e+01 -2.022 0.0545 .
## ImageSetopioid -7.434e-02 1.518e-01 6.562e+03 -0.490 0.6243
## Time:VisitV2 -4.095e-03 4.632e-04 6.562e+03 -8.842 < 2e-16 ***
## Age:SexMale 1.661e-01 6.960e-02 2.410e+01 2.386 0.0252 *
## Age:ImageSetopioid 2.162e-03 3.915e-03 6.562e+03 0.552 0.5809
## SexMale:ImageSetopioid 5.365e-01 2.421e-01 6.562e+03 2.216 0.0267 *
## Age:SexMale:ImageSetopioid -1.445e-02 6.273e-03 6.562e+03 -2.303 0.0213 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Time VistV2 Age SexMal ImgStp Tm:VV2 Ag:SxM Ag:ImS SxM:IS
## Time -0.018
## VisitV2 -0.013 0.602
## Age -0.969 0.000 -0.001
## SexMale -0.629 0.000 0.000 0.610
## ImageSetopd -0.045 0.001 -0.013 0.044 0.028
## Time:VistV2 0.011 -0.697 -0.865 0.001 0.001 0.015
## Age:SexMale 0.606 0.000 0.000 -0.625 -0.979 -0.027 -0.001
## Age:ImgStpd 0.044 0.002 0.018 -0.045 -0.027 -0.969 -0.021 0.028
## SxMl:ImgStp 0.028 0.005 0.012 -0.027 -0.045 -0.627 -0.013 0.044 0.608
## Ag:SxMl:ImS -0.027 -0.006 -0.016 0.028 0.044 0.605 0.018 -0.045 -0.624 -0.979
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 0.006 0.006 1 6562.0 0.0058 0.93939
## Visit 62.493 62.493 1 6562.3 62.8088 2.657e-15 ***
## Age 4.898 4.898 1 24.0 4.9228 0.03621 *
## Sex 3.682 3.682 1 24.0 3.7001 0.06635 .
## ImageSet 2.553 2.553 1 6562.0 2.5660 0.10923
## Time:Visit 77.782 77.782 1 6562.0 78.1749 < 2.2e-16 ***
## Age:Sex 5.193 5.193 1 24.0 5.2193 0.03147 *
## Age:ImageSet 2.593 2.593 1 6562.0 2.6062 0.10649
## Sex:ImageSet 4.886 4.886 1 6562.0 4.9106 0.02673 *
## Age:Sex:ImageSet 5.279 5.279 1 6562.0 5.3056 0.02129 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'arousal')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 20211.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.0751 -0.3333 0.0039 0.4164 7.0093
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 4.586 2.141
## Residual 1.206 1.098
## Number of obs: 6600, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.884e+00 2.526e+00 2.406e+01 1.538 0.13714
## Time -6.623e-04 3.556e-04 6.565e+03 -1.863 0.06256 .
## VisitV2 -1.403e-01 5.441e-02 6.565e+03 -2.578 0.00995 **
## Age -1.766e-02 6.515e-02 2.405e+01 -0.271 0.78866
## SexMale -2.643e+00 4.016e+00 2.405e+01 -0.658 0.51677
## ImageSetopioid 1.941e-02 1.671e-01 6.565e+03 0.116 0.90751
## Time:VisitV2 8.060e-04 5.098e-04 6.565e+03 1.581 0.11393
## Age:SexMale 1.414e-01 1.042e-01 2.405e+01 1.357 0.18747
## Age:ImageSetopioid 1.039e-03 4.310e-03 6.565e+03 0.241 0.80946
## SexMale:ImageSetopioid 3.256e-01 2.665e-01 6.565e+03 1.222 0.22190
## Age:SexMale:ImageSetopioid -1.157e-02 6.906e-03 6.565e+03 -1.675 0.09400 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Time VistV2 Age SexMal ImgStp Tm:VV2 Ag:SxM Ag:ImS SxM:IS
## Time -0.013
## VisitV2 -0.010 0.602
## Age -0.969 0.000 -0.001
## SexMale -0.629 0.000 0.000 0.610
## ImageSetopd -0.033 0.000 -0.013 0.032 0.021
## Time:VistV2 0.008 -0.698 -0.865 0.001 0.000 0.015
## Age:SexMale 0.606 0.000 0.000 -0.625 -0.979 -0.020 -0.001
## Age:ImgStpd 0.032 0.003 0.019 -0.033 -0.020 -0.970 -0.022 0.021
## SxMl:ImgStp 0.021 0.005 0.012 -0.020 -0.033 -0.627 -0.014 0.032 0.608
## Ag:SxMl:ImS -0.020 -0.006 -0.016 0.021 0.032 0.605 0.018 -0.033 -0.624 -0.979
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 1.2487 1.2487 1 6565.0 1.0354 0.308928
## Visit 8.0167 8.0167 1 6565.2 6.6471 0.009953 **
## Age 1.1415 1.1415 1 24.0 0.9465 0.340326
## Sex 0.4604 0.4604 1 24.0 0.3817 0.542509
## ImageSet 2.2546 2.2546 1 6565.0 1.8694 0.171590
## Time:Visit 3.0144 3.0144 1 6565.0 2.4994 0.113935
## Age:Sex 2.0443 2.0443 1 24.0 1.6950 0.205298
## Age:ImageSet 2.2771 2.2771 1 6565.0 1.8880 0.169469
## Sex:ImageSet 1.7998 1.7998 1 6565.0 1.4923 0.221897
## Age:Sex:ImageSet 3.3833 3.3833 1 6565.0 2.8053 0.094003 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'typicality')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 52301.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9362 -0.1409 0.1529 0.4851 3.6520
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 139 11.79
## Residual 158 12.57
## Number of obs: 6600, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.726e+01 1.396e+01 2.427e+01 5.533 1.04e-05 ***
## Time -2.308e-02 4.071e-03 6.565e+03 -5.670 1.49e-08 ***
## VisitV2 -2.251e+00 6.229e-01 6.566e+03 -3.614 0.000304 ***
## Age 4.869e-01 3.600e-01 2.423e+01 1.352 0.188721
## SexMale 3.476e+00 2.219e+01 2.423e+01 0.157 0.876857
## ImageSetopioid 6.289e-01 1.912e+00 6.565e+03 0.329 0.742270
## Time:VisitV2 -2.672e-02 5.836e-03 6.565e+03 -4.579 4.77e-06 ***
## Age:SexMale -1.764e-01 5.760e-01 2.423e+01 -0.306 0.761983
## Age:ImageSetopioid 4.352e-03 4.934e-02 6.565e+03 0.088 0.929714
## SexMale:ImageSetopioid -2.281e+00 3.051e+00 6.565e+03 -0.748 0.454759
## Age:SexMale:ImageSetopioid 3.339e-02 7.906e-02 6.565e+03 0.422 0.672797
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Time VistV2 Age SexMal ImgStp Tm:VV2 Ag:SxM Ag:ImS SxM:IS
## Time -0.027
## VisitV2 -0.020 0.602
## Age -0.969 0.000 -0.002
## SexMale -0.629 0.000 0.000 0.610
## ImageSetopd -0.069 0.000 -0.013 0.066 0.043
## Time:VistV2 0.017 -0.698 -0.865 0.002 0.001 0.015
## Age:SexMale 0.606 0.000 0.001 -0.625 -0.979 -0.042 -0.001
## Age:ImgStpd 0.066 0.003 0.019 -0.069 -0.042 -0.970 -0.022 0.043
## SxMl:ImgStp 0.043 0.005 0.012 -0.042 -0.069 -0.627 -0.014 0.067 0.608
## Ag:SxMl:ImS -0.041 -0.006 -0.016 0.043 0.067 0.605 0.018 -0.069 -0.624 -0.979
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 24654.7 24654.7 1 6565.1 155.9937 < 2.2e-16 ***
## Visit 2064.0 2064.0 1 6565.8 13.0590 0.0003041 ***
## Age 320.6 320.6 1 24.0 2.0286 0.1672316
## Sex 1.8 1.8 1 24.0 0.0111 0.9168811
## ImageSet 17.8 17.8 1 6565.0 0.1124 0.7374118
## Time:Visit 3313.3 3313.3 1 6565.1 20.9635 4.767e-06 ***
## Age:Sex 12.2 12.2 1 24.0 0.0773 0.7834056
## Age:ImageSet 44.8 44.8 1 6565.0 0.2835 0.5944163
## Sex:ImageSet 88.3 88.3 1 6565.0 0.5588 0.4547588
## Age:Sex:ImageSet 28.2 28.2 1 6565.0 0.1784 0.6727975
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'craving')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 53717
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.2791 -0.1078 0.1569 0.4930 3.2842
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 130.9 11.44
## Residual 196.1 14.00
## Number of obs: 6600, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.594e+01 1.357e+01 2.435e+01 5.596 8.82e-06 ***
## Time -2.846e-02 4.535e-03 6.565e+03 -6.276 3.69e-10 ***
## VisitV2 -1.413e+00 6.939e-01 6.566e+03 -2.037 0.0417 *
## Age 4.763e-01 3.499e-01 2.430e+01 1.361 0.1859
## SexMale 1.086e+00 2.157e+01 2.430e+01 0.050 0.9603
## ImageSetopioid 5.836e-01 2.130e+00 6.565e+03 0.274 0.7841
## Time:VisitV2 -2.851e-02 6.501e-03 6.565e+03 -4.386 1.17e-05 ***
## Age:SexMale -1.001e-01 5.597e-01 2.430e+01 -0.179 0.8596
## Age:ImageSetopioid 1.003e-02 5.496e-02 6.565e+03 0.183 0.8552
## SexMale:ImageSetopioid -3.359e+00 3.399e+00 6.565e+03 -0.988 0.3231
## Age:SexMale:ImageSetopioid 5.819e-02 8.807e-02 6.565e+03 0.661 0.5088
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Time VistV2 Age SexMal ImgStp Tm:VV2 Ag:SxM Ag:ImS SxM:IS
## Time -0.031
## VisitV2 -0.023 0.602
## Age -0.969 0.000 -0.002
## SexMale -0.628 0.000 0.000 0.610
## ImageSetopd -0.079 0.000 -0.013 0.076 0.049
## Time:VistV2 0.019 -0.698 -0.865 0.003 0.001 0.015
## Age:SexMale 0.606 0.000 0.001 -0.625 -0.979 -0.048 -0.001
## Age:ImgStpd 0.076 0.003 0.019 -0.079 -0.048 -0.970 -0.022 0.049
## SxMl:ImgStp 0.049 0.005 0.012 -0.048 -0.079 -0.627 -0.014 0.077 0.608
## Ag:SxMl:ImS -0.047 -0.006 -0.016 0.049 0.077 0.605 0.018 -0.079 -0.624 -0.979
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 33878 33878 1 6565.1 172.7388 < 2.2e-16 ***
## Visit 814 814 1 6566.0 4.1480 0.04172 *
## Age 501 501 1 24.0 2.5536 0.12313
## Sex 0 0 1 24.0 0.0008 0.97822
## ImageSet 82 82 1 6565.0 0.4158 0.51909
## Time:Visit 3773 3773 1 6565.1 19.2366 1.173e-05 ***
## Age:Sex 3 3 1 24.0 0.0162 0.89985
## Age:ImageSet 155 155 1 6565.0 0.7897 0.37423
## Sex:ImageSet 192 192 1 6565.0 0.9765 0.32309
## Age:Sex:ImageSet 86 86 1 6565.0 0.4366 0.50879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_neutral(subject_data[subject_data$ImageSet %in% c('control'),], 'valence')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex) + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 9141.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2538 -0.1570 0.0124 0.0963 6.5360
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.0948 1.0463
## Residual 0.8862 0.9414
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.347e+00 1.237e+00 2.406e+01 3.513 0.00178 **
## Time 1.159e-03 4.728e-04 3.269e+03 2.452 0.01426 *
## VisitV2 3.257e-02 6.516e-02 3.269e+03 0.500 0.61719
## Age 1.340e-02 3.190e-02 2.400e+01 0.420 0.67822
## SexMale 1.661e+00 1.967e+00 2.400e+01 0.845 0.40664
## Time:VisitV2 -8.366e-04 6.698e-04 3.269e+03 -1.249 0.21176
## Age:SexMale -4.637e-02 5.104e-02 2.400e+01 -0.909 0.37264
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Time VistV2 Age SexMal Tm:VV2
## Time -0.032
## VisitV2 -0.028 0.612
## Age -0.969 0.000 0.002
## SexMale -0.628 -0.001 0.001 0.610
## Time:VistV2 0.025 -0.706 -0.862 -0.002 0.000
## Age:SexMale 0.606 0.000 -0.001 -0.625 -0.979 0.000
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 4.3331 4.3331 1 3269.4 4.8898 0.02709 *
## Visit 0.2214 0.2214 1 3269.4 0.2499 0.61719
## Age 0.1303 0.1303 1 24.0 0.1470 0.70475
## Sex 0.6322 0.6322 1 24.0 0.7135 0.40664
## Time:Visit 1.3824 1.3824 1 3269.3 1.5600 0.21176
## Age:Sex 0.7314 0.7314 1 24.0 0.8254 0.37264
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_neutral(subject_data[subject_data$ImageSet %in% c('control'),], 'arousal')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex) + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 11260.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9989 -0.3580 -0.0204 0.2142 5.9661
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.828 1.352
## Residual 1.688 1.299
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.039e+00 1.600e+00 2.403e+01 1.275 0.21467
## Time -8.255e-05 6.525e-04 3.269e+03 -0.127 0.89933
## VisitV2 -2.582e-02 8.993e-02 3.269e+03 -0.287 0.77405
## Age -1.185e-02 4.124e-02 2.396e+01 -0.287 0.77633
## SexMale 2.685e+00 2.543e+00 2.397e+01 1.056 0.30153
## Time:VisitV2 2.666e-03 9.244e-04 3.269e+03 2.884 0.00395 **
## Age:SexMale -3.453e-02 6.598e-02 2.397e+01 -0.523 0.60554
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Time VistV2 Age SexMal Tm:VV2
## Time -0.034
## VisitV2 -0.030 0.612
## Age -0.969 0.000 0.002
## SexMale -0.628 -0.001 0.001 0.610
## Time:VistV2 0.026 -0.706 -0.862 -0.002 0.000
## Age:SexMale 0.606 0.001 -0.001 -0.625 -0.979 0.000
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 12.3426 12.3426 1 3269.4 7.3127 0.006882 **
## Visit 0.1391 0.1391 1 3269.4 0.0824 0.774055
## Age 1.3146 1.3146 1 24.0 0.7789 0.386259
## Sex 1.8819 1.8819 1 24.0 1.1150 0.301526
## Time:Visit 14.0390 14.0390 1 3269.3 8.3178 0.003951 **
## Age:Sex 0.4623 0.4623 1 24.0 0.2739 0.605545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_neutral(subject_data[subject_data$ImageSet %in% c('control'),], 'typicality')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex) + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 27913.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5429 -0.3959 -0.1737 0.0319 6.0613
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 179.3 13.39
## Residual 266.2 16.31
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.569e+01 1.589e+01 2.409e+01 2.247 0.034098 *
## Time -1.387e-02 8.194e-03 3.270e+03 -1.693 0.090569 .
## VisitV2 -4.209e+00 1.129e+00 3.270e+03 -3.727 0.000197 ***
## Age -6.669e-01 4.094e-01 2.399e+01 -1.629 0.116328
## SexMale -1.732e+01 2.524e+01 2.399e+01 -0.686 0.499188
## Time:VisitV2 5.767e-02 1.161e-02 3.270e+03 4.968 7.12e-07 ***
## Age:SexMale 5.726e-01 6.549e-01 2.399e+01 0.874 0.390600
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Time VistV2 Age SexMal Tm:VV2
## Time -0.043
## VisitV2 -0.038 0.612
## Age -0.968 0.000 0.002
## SexMale -0.628 -0.001 0.001 0.610
## Time:VistV2 0.033 -0.706 -0.862 -0.003 0.000
## Age:SexMale 0.605 0.001 -0.001 -0.625 -0.979 0.000
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 1767.2 1767.2 1 3269.6 6.6393 0.0100189 *
## Visit 3697.0 3697.0 1 3269.6 13.8894 0.0001972 ***
## Age 359.6 359.6 1 24.0 1.3510 0.2565287
## Sex 125.3 125.3 1 24.0 0.4708 0.4991875
## Time:Visit 6568.8 6568.8 1 3269.6 24.6784 7.122e-07 ***
## Age:Sex 203.5 203.5 1 24.0 0.7645 0.3905997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_neutral(subject_data[subject_data$ImageSet %in% c('control'),], 'craving')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex) + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 28253.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6918 -0.3422 -0.1704 -0.0075 5.7139
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 197.7 14.06
## Residual 295.1 17.18
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.682e+01 1.669e+01 2.409e+01 1.607 0.12102
## Time -9.995e-03 8.627e-03 3.270e+03 -1.159 0.24670
## VisitV2 -2.631e+00 1.189e+00 3.270e+03 -2.212 0.02700 *
## Age -4.613e-01 4.299e-01 2.399e+01 -1.073 0.29397
## SexMale -1.198e+01 2.651e+01 2.399e+01 -0.452 0.65528
## Time:VisitV2 3.454e-02 1.222e-02 3.270e+03 2.826 0.00475 **
## Age:SexMale 4.437e-01 6.879e-01 2.399e+01 0.645 0.52506
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Time VistV2 Age SexMal Tm:VV2
## Time -0.043
## VisitV2 -0.038 0.612
## Age -0.968 0.000 0.002
## SexMale -0.628 -0.001 0.001 0.610
## Time:VistV2 0.033 -0.706 -0.862 -0.003 0.000
## Age:SexMale 0.605 0.001 -0.001 -0.625 -0.979 0.000
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 417.58 417.58 1 3269.6 1.4151 0.234296
## Visit 1444.44 1444.44 1 3269.6 4.8950 0.027004 *
## Age 143.07 143.07 1 24.0 0.4848 0.492929
## Sex 60.30 60.30 1 24.0 0.2044 0.655284
## Time:Visit 2356.19 2356.19 1 3269.6 7.9847 0.004746 **
## Age:Sex 122.76 122.76 1 24.0 0.4160 0.525059
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category <- function(dset, RatingType){
this_lme <- lmer(Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id),
data = dset[dset$RatingType == RatingType,])
print(plot_model(this_lme, title = paste(dset$ImageSet[1], RatingType), show.values = TRUE, type = 'std'))
print(summary(this_lme))
#for using nlme instead of lme4--but can't get standardized betas in the forest plot
#that_lme <- lme(fixed = Value ~ Time * Visit + Age + Sex + ImageSet, random = ~ 1 | id,
# data = dset[dset$RatingType == RatingType,])
#print(summary(that_lme))
print(anova(this_lme))
}
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('meth'),], 'valence')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 9492.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8790 -0.4014 0.0437 0.3901 6.7604
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.7891 1.3376
## Residual 0.9795 0.9897
## Number of obs: 3297, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.274e+00 1.599e+00 2.314e+01 2.048 0.05212
## Time 1.955e-03 4.525e-04 3.256e+03 4.320 1.61e-05
## VisitV2 4.048e-01 6.986e-02 3.256e+03 5.793 7.55e-09
## Age -1.740e-02 4.121e-02 2.300e+01 -0.422 0.67668
## SexMale -7.375e+00 2.657e+00 2.300e+01 -2.775 0.01076
## InjectionUserTRUE 1.532e+00 6.288e-01 2.394e+01 2.436 0.02268
## Categorymeth_face_activities 6.199e-02 1.183e-01 3.256e+03 0.524 0.60033
## Categorymeth_hand 2.668e-01 1.366e-01 3.256e+03 1.953 0.05087
## Categorymeth_injection_hand -1.199e-01 1.366e-01 3.256e+03 -0.878 0.38014
## Categorymeth_instrument 1.384e-01 1.044e-01 3.256e+03 1.325 0.18518
## Categorymeth_instrument_hand 6.609e-03 1.183e-01 3.256e+03 0.056 0.95546
## Time:VisitV2 -4.202e-03 6.526e-04 3.256e+03 -6.439 1.38e-10
## Age:SexMale 2.179e-01 6.911e-02 2.300e+01 3.153 0.00445
## InjectionUserTRUE:Categorymeth_face_activities -6.608e-02 1.370e-01 3.256e+03 -0.482 0.62971
## InjectionUserTRUE:Categorymeth_hand -2.483e-01 1.583e-01 3.256e+03 -1.569 0.11684
## InjectionUserTRUE:Categorymeth_injection_hand -2.217e-01 1.583e-01 3.256e+03 -1.401 0.16129
## InjectionUserTRUE:Categorymeth_instrument -1.966e-01 1.209e-01 3.256e+03 -1.626 0.10408
## InjectionUserTRUE:Categorymeth_instrument_hand -3.380e-02 1.370e-01 3.256e+03 -0.247 0.80517
##
## (Intercept) .
## Time ***
## VisitV2 ***
## Age
## SexMale *
## InjectionUserTRUE *
## Categorymeth_face_activities
## Categorymeth_hand .
## Categorymeth_injection_hand
## Categorymeth_instrument
## Categorymeth_instrument_hand
## Time:VisitV2 ***
## Age:SexMale **
## InjectionUserTRUE:Categorymeth_face_activities
## InjectionUserTRUE:Categorymeth_hand
## InjectionUserTRUE:Categorymeth_injection_hand
## InjectionUserTRUE:Categorymeth_instrument
## InjectionUserTRUE:Categorymeth_instrument_hand
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 0.197 0.197 1 3256.2 0.2008 0.654106
## Visit 32.875 32.875 1 3256.4 33.5641 7.552e-09 ***
## Age 7.529 7.529 1 23.0 7.6871 0.010831 *
## Sex 7.545 7.545 1 23.0 7.7028 0.010760 *
## InjectionUser 4.978 4.978 1 23.0 5.0823 0.033986 *
## Category 19.930 3.986 5 3256.0 4.0695 0.001101 **
## Time:Visit 40.613 40.613 1 3256.1 41.4637 1.377e-10 ***
## Age:Sex 9.736 9.736 1 23.0 9.9405 0.004452 **
## InjectionUser:Category 5.352 1.070 5 3256.0 1.0927 0.362171
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('meth'),], 'arousal')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 10149.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3380 -0.3677 0.0272 0.3884 6.8963
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 4.323 2.079
## Residual 1.188 1.090
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.220e+00 2.481e+00 2.307e+01 1.298 0.2072
## Time -2.886e-04 4.981e-04 3.259e+03 -0.579 0.5623
## VisitV2 -7.477e-02 7.691e-02 3.259e+03 -0.972 0.3310
## Age -3.242e-02 6.398e-02 2.300e+01 -0.507 0.6172
## SexMale -4.722e+00 4.126e+00 2.300e+01 -1.145 0.2642
## InjectionUserTRUE 1.637e+00 9.716e-01 2.347e+01 1.685 0.1052
## Categorymeth_face_activities 1.790e-01 1.303e-01 3.259e+03 1.374 0.1694
## Categorymeth_hand 1.521e-01 1.504e-01 3.259e+03 1.011 0.3121
## Categorymeth_injection_hand -1.696e-02 1.504e-01 3.259e+03 -0.113 0.9102
## Categorymeth_instrument -3.046e-02 1.149e-01 3.259e+03 -0.265 0.7910
## Categorymeth_instrument_hand 7.707e-02 1.303e-01 3.259e+03 0.592 0.5542
## Time:VisitV2 5.338e-04 7.183e-04 3.259e+03 0.743 0.4574
## Age:SexMale 1.965e-01 1.073e-01 2.300e+01 1.831 0.0801 .
## InjectionUserTRUE:Categorymeth_face_activities -3.047e-01 1.509e-01 3.259e+03 -2.020 0.0435 *
## InjectionUserTRUE:Categorymeth_hand -2.955e-01 1.743e-01 3.259e+03 -1.696 0.0900 .
## InjectionUserTRUE:Categorymeth_injection_hand 1.375e-01 1.743e-01 3.259e+03 0.789 0.4303
## InjectionUserTRUE:Categorymeth_instrument -1.439e-01 1.331e-01 3.259e+03 -1.081 0.2797
## InjectionUserTRUE:Categorymeth_instrument_hand -2.292e-01 1.509e-01 3.259e+03 -1.519 0.1289
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 0.0043 0.0043 1 3259.1 0.0037 0.95182
## Visit 1.1225 1.1225 1 3259.2 0.9451 0.33103
## Age 1.9568 1.9568 1 23.0 1.6476 0.21206
## Sex 1.5558 1.5558 1 23.0 1.3100 0.26417
## InjectionUser 2.8509 2.8509 1 23.0 2.4005 0.13494
## Category 7.6126 1.5225 5 3259.0 1.2820 0.26865
## Time:Visit 0.6560 0.6560 1 3259.1 0.5524 0.45741
## Age:Sex 3.9806 3.9806 1 23.0 3.3517 0.08012 .
## InjectionUser:Category 12.0005 2.4001 5 3259.0 2.0209 0.07263 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('meth'),], 'typicality')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 26292.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.8094 -0.1627 0.1535 0.4861 2.1803
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 152.2 12.33
## Residual 164.0 12.81
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.795e+01 1.480e+01 2.328e+01 5.267 2.33e-05
## Time -2.021e-02 5.854e-03 3.259e+03 -3.453 0.000562
## VisitV2 -2.037e+00 9.038e-01 3.260e+03 -2.254 0.024275
## Age 5.097e-01 3.808e-01 2.300e+01 1.338 0.193879
## SexMale 6.700e+00 2.456e+01 2.300e+01 0.273 0.787409
## InjectionUserTRUE -8.166e-01 5.867e+00 2.485e+01 -0.139 0.890414
## Categorymeth_face_activities 1.212e+00 1.531e+00 3.259e+03 0.792 0.428425
## Categorymeth_hand 1.153e+00 1.768e+00 3.259e+03 0.652 0.514150
## Categorymeth_injection_hand 1.782e+00 1.768e+00 3.259e+03 1.008 0.313462
## Categorymeth_instrument -1.790e+00 1.351e+00 3.259e+03 -1.325 0.185166
## Categorymeth_instrument_hand 7.153e-02 1.531e+00 3.259e+03 0.047 0.962741
## Time:VisitV2 -3.023e-02 8.441e-03 3.259e+03 -3.581 0.000348
## Age:SexMale -2.617e-01 6.387e-01 2.300e+01 -0.410 0.685802
## InjectionUserTRUE:Categorymeth_face_activities -1.685e+00 1.773e+00 3.259e+03 -0.950 0.342153
## InjectionUserTRUE:Categorymeth_hand -2.984e+00 2.048e+00 3.259e+03 -1.457 0.145144
## InjectionUserTRUE:Categorymeth_injection_hand -3.030e+00 2.048e+00 3.259e+03 -1.479 0.139136
## InjectionUserTRUE:Categorymeth_instrument -4.023e-01 1.564e+00 3.259e+03 -0.257 0.796991
## InjectionUserTRUE:Categorymeth_instrument_hand -2.953e+00 1.773e+00 3.259e+03 -1.665 0.095947
##
## (Intercept) ***
## Time ***
## VisitV2 *
## Age
## SexMale
## InjectionUserTRUE
## Categorymeth_face_activities
## Categorymeth_hand
## Categorymeth_injection_hand
## Categorymeth_instrument
## Categorymeth_instrument_hand
## Time:VisitV2 ***
## Age:SexMale
## InjectionUserTRUE:Categorymeth_face_activities
## InjectionUserTRUE:Categorymeth_hand
## InjectionUserTRUE:Categorymeth_injection_hand
## InjectionUserTRUE:Categorymeth_instrument
## InjectionUserTRUE:Categorymeth_instrument_hand .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 11466.4 11466.4 1 3259.3 69.9058 < 2.2e-16 ***
## Visit 833.2 833.2 1 3259.8 5.0796 0.0242748 *
## Age 252.8 252.8 1 23.0 1.5409 0.2269949
## Sex 12.2 12.2 1 23.0 0.0744 0.7874094
## InjectionUser 35.0 35.0 1 23.1 0.2133 0.6485393
## Category 2389.7 477.9 5 3259.0 2.9138 0.0125026 *
## Time:Visit 2103.2 2103.2 1 3259.2 12.8223 0.0003475 ***
## Age:Sex 27.5 27.5 1 23.0 0.1679 0.6858018
## InjectionUser:Category 968.8 193.8 5 3259.0 1.1813 0.3156999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('meth'),], 'craving')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 26970.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.9852 -0.1233 0.1477 0.5068 3.0025
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 140.9 11.87
## Residual 202.1 14.22
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.650e+01 1.428e+01 2.337e+01 5.358 1.83e-05
## Time -2.095e-02 6.498e-03 3.259e+03 -3.225 0.00127
## VisitV2 1.267e-01 1.003e+00 3.260e+03 0.126 0.89953
## Age 4.821e-01 3.670e-01 2.300e+01 1.313 0.20198
## SexMale 2.004e+00 2.367e+01 2.300e+01 0.085 0.93326
## InjectionUserTRUE 2.737e-01 5.689e+00 2.548e+01 0.048 0.96200
## Categorymeth_face_activities 2.876e-01 1.699e+00 3.259e+03 0.169 0.86561
## Categorymeth_hand 2.905e-02 1.962e+00 3.259e+03 0.015 0.98819
## Categorymeth_injection_hand 2.994e-01 1.962e+00 3.259e+03 0.153 0.87874
## Categorymeth_instrument -3.114e+00 1.499e+00 3.259e+03 -2.077 0.03784
## Categorymeth_instrument_hand -1.367e+00 1.700e+00 3.259e+03 -0.805 0.42115
## Time:VisitV2 -4.183e-02 9.370e-03 3.259e+03 -4.465 8.29e-06
## Age:SexMale -1.245e-01 6.156e-01 2.300e+01 -0.202 0.84147
## InjectionUserTRUE:Categorymeth_face_activities -2.100e+00 1.968e+00 3.259e+03 -1.067 0.28604
## InjectionUserTRUE:Categorymeth_hand -2.694e+00 2.273e+00 3.259e+03 -1.185 0.23604
## InjectionUserTRUE:Categorymeth_injection_hand -1.215e+00 2.273e+00 3.259e+03 -0.534 0.59318
## InjectionUserTRUE:Categorymeth_instrument 2.638e-01 1.736e+00 3.259e+03 0.152 0.87923
## InjectionUserTRUE:Categorymeth_instrument_hand -2.108e+00 1.968e+00 3.259e+03 -1.071 0.28414
##
## (Intercept) ***
## Time **
## VisitV2
## Age
## SexMale
## InjectionUserTRUE
## Categorymeth_face_activities
## Categorymeth_hand
## Categorymeth_injection_hand
## Categorymeth_instrument *
## Categorymeth_instrument_hand
## Time:VisitV2 ***
## Age:SexMale
## InjectionUserTRUE:Categorymeth_face_activities
## InjectionUserTRUE:Categorymeth_hand
## InjectionUserTRUE:Categorymeth_injection_hand
## InjectionUserTRUE:Categorymeth_instrument
## InjectionUserTRUE:Categorymeth_instrument_hand
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 16109.0 16109.0 1 3259.4 79.7049 < 2.2e-16 ***
## Visit 3.2 3.2 1 3260.1 0.0159 0.899528
## Age 411.8 411.8 1 23.0 2.0375 0.166898
## Sex 1.4 1.4 1 23.0 0.0072 0.933264
## InjectionUser 7.0 7.0 1 23.1 0.0348 0.853668
## Category 3508.8 701.8 5 3259.0 3.4722 0.003926 **
## Time:Visit 4028.7 4028.7 1 3259.3 19.9333 8.291e-06 ***
## Age:Sex 8.3 8.3 1 23.0 0.0409 0.841467
## InjectionUser:Category 840.6 168.1 5 3259.0 0.8318 0.526839
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('opioid'),], 'valence')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 9405.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0560 -0.3698 0.0620 0.4468 4.4949
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.8191 1.3487
## Residual 0.9511 0.9753
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.471e+00 1.612e+00 2.313e+01 2.153
## Time 2.090e-03 4.501e-04 3.259e+03 4.642
## VisitV2 3.703e-01 6.814e-02 3.259e+03 5.435
## Age -1.310e-02 4.154e-02 2.300e+01 -0.315
## SexMale -6.565e+00 2.679e+00 2.300e+01 -2.451
## InjectionUserTRUE 1.321e+00 6.337e-01 2.390e+01 2.084
## Categoryopioid_face_activities -2.766e-01 1.166e-01 3.259e+03 -2.373
## Categoryopioid_hand 2.285e-01 1.346e-01 3.259e+03 1.697
## Categoryopioid_injection_hand -3.905e-01 1.346e-01 3.259e+03 -2.901
## Categoryopioid_injection_instrument -1.485e-01 1.028e-01 3.259e+03 -1.444
## Categoryopioid_instrument_hand -1.613e-01 1.166e-01 3.259e+03 -1.384
## Time:VisitV2 -3.894e-03 6.415e-04 3.259e+03 -6.070
## Age:SexMale 1.960e-01 6.968e-02 2.300e+01 2.814
## InjectionUserTRUE:Categoryopioid_face_activities -1.016e-01 1.350e-01 3.259e+03 -0.753
## InjectionUserTRUE:Categoryopioid_hand -2.469e-01 1.559e-01 3.259e+03 -1.584
## InjectionUserTRUE:Categoryopioid_injection_hand -1.306e-01 1.559e-01 3.259e+03 -0.838
## InjectionUserTRUE:Categoryopioid_injection_instrument -1.424e-01 1.191e-01 3.259e+03 -1.196
## InjectionUserTRUE:Categoryopioid_instrument_hand -1.607e-01 1.351e-01 3.259e+03 -1.190
## Pr(>|t|)
## (Intercept) 0.04200 *
## Time 3.58e-06 ***
## VisitV2 5.89e-08 ***
## Age 0.75535
## SexMale 0.02228 *
## InjectionUserTRUE 0.04803 *
## Categoryopioid_face_activities 0.01771 *
## Categoryopioid_hand 0.08975 .
## Categoryopioid_injection_hand 0.00374 **
## Categoryopioid_injection_instrument 0.14885
## Categoryopioid_instrument_hand 0.16660
## Time:VisitV2 1.43e-09 ***
## Age:SexMale 0.00986 **
## InjectionUserTRUE:Categoryopioid_face_activities 0.45173
## InjectionUserTRUE:Categoryopioid_hand 0.11336
## InjectionUserTRUE:Categoryopioid_injection_hand 0.40224
## InjectionUserTRUE:Categoryopioid_injection_instrument 0.23178
## InjectionUserTRUE:Categoryopioid_instrument_hand 0.23413
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 0.188 0.188 1 3259.1 0.1982 0.656224
## Visit 28.093 28.093 1 3259.4 29.5366 5.891e-08 ***
## Age 6.189 6.189 1 23.0 6.5076 0.017859 *
## Sex 5.712 5.712 1 23.0 6.0059 0.022281 *
## InjectionUser 3.418 3.418 1 23.0 3.5939 0.070611 .
## Category 64.145 12.829 5 3259.0 13.4883 4.779e-13 ***
## Time:Visit 35.042 35.042 1 3259.2 36.8430 1.428e-09 ***
## Age:Sex 7.529 7.529 1 23.0 7.9159 0.009859 **
## InjectionUser:Category 2.809 0.562 5 3259.0 0.5907 0.707114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('opioid'),], 'arousal')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 10093.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.4731 -0.2912 0.0237 0.4270 3.5615
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 4.272 2.067
## Residual 1.167 1.080
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.271e+00 2.466e+00 2.307e+01 1.326
## Time -1.274e-03 4.987e-04 3.259e+03 -2.555
## VisitV2 -2.422e-01 7.550e-02 3.259e+03 -3.208
## Age -3.333e-02 6.360e-02 2.300e+01 -0.524
## SexMale -4.659e+00 4.101e+00 2.300e+01 -1.136
## InjectionUserTRUE 1.557e+00 9.658e-01 2.347e+01 1.612
## Categoryopioid_face_activities 1.214e-02 1.292e-01 3.259e+03 0.094
## Categoryopioid_hand 1.236e-01 1.491e-01 3.259e+03 0.829
## Categoryopioid_injection_hand 3.257e-02 1.491e-01 3.259e+03 0.218
## Categoryopioid_injection_instrument 7.458e-02 1.139e-01 3.259e+03 0.655
## Categoryopioid_instrument_hand 2.097e-01 1.292e-01 3.259e+03 1.623
## Time:VisitV2 1.486e-03 7.107e-04 3.259e+03 2.092
## Age:SexMale 1.919e-01 1.067e-01 2.300e+01 1.799
## InjectionUserTRUE:Categoryopioid_face_activities 2.936e-01 1.496e-01 3.259e+03 1.963
## InjectionUserTRUE:Categoryopioid_hand -1.746e-01 1.727e-01 3.259e+03 -1.011
## InjectionUserTRUE:Categoryopioid_injection_hand 4.901e-01 1.727e-01 3.259e+03 2.838
## InjectionUserTRUE:Categoryopioid_injection_instrument 1.494e-01 1.319e-01 3.259e+03 1.132
## InjectionUserTRUE:Categoryopioid_instrument_hand -1.336e-03 1.496e-01 3.259e+03 -0.009
## Pr(>|t|)
## (Intercept) 0.19775
## Time 0.01066 *
## VisitV2 0.00135 **
## Age 0.60527
## SexMale 0.26765
## InjectionUserTRUE 0.12039
## Categoryopioid_face_activities 0.92511
## Categoryopioid_hand 0.40723
## Categoryopioid_injection_hand 0.82711
## Categoryopioid_injection_instrument 0.51274
## Categoryopioid_instrument_hand 0.10474
## Time:VisitV2 0.03656 *
## Age:SexMale 0.08515 .
## InjectionUserTRUE:Categoryopioid_face_activities 0.04977 *
## InjectionUserTRUE:Categoryopioid_hand 0.31206
## InjectionUserTRUE:Categoryopioid_injection_hand 0.00457 **
## InjectionUserTRUE:Categoryopioid_injection_instrument 0.25759
## InjectionUserTRUE:Categoryopioid_instrument_hand 0.99288
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 2.6064 2.6064 1 3259.1 2.2326 0.135222
## Visit 12.0177 12.0177 1 3259.2 10.2943 0.001347 **
## Age 1.7628 1.7628 1 23.0 1.5100 0.231558
## Sex 1.5066 1.5066 1 23.0 1.2905 0.267647
## InjectionUser 3.5789 3.5789 1 23.0 3.0657 0.093281 .
## Category 17.8034 3.5607 5 3259.0 3.0501 0.009456 **
## Time:Visit 5.1069 5.1069 1 3259.1 4.3745 0.036557 *
## Age:Sex 3.7784 3.7784 1 23.0 3.2365 0.085149 .
## InjectionUser:Category 20.3376 4.0675 5 3259.0 3.4842 0.003828 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('opioid'),], 'typicality')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 25997.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.8348 -0.1406 0.1551 0.4661 3.8109
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 136.3 11.67
## Residual 150.0 12.25
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 7.898e+01 1.401e+01 2.328e+01 5.638
## Time -2.598e-02 5.652e-03 3.259e+03 -4.596
## VisitV2 -2.593e+00 8.556e-01 3.260e+03 -3.030
## Age 5.158e-01 3.604e-01 2.300e+01 1.431
## SexMale 4.633e+00 2.324e+01 2.300e+01 0.199
## InjectionUserTRUE -4.029e+00 5.555e+00 2.489e+01 -0.725
## Categoryopioid_face_activities 1.616e+00 1.464e+00 3.259e+03 1.104
## Categoryopioid_hand -1.970e+00 1.690e+00 3.259e+03 -1.165
## Categoryopioid_injection_hand 1.826e+00 1.690e+00 3.259e+03 1.080
## Categoryopioid_injection_instrument -1.261e+00 1.291e+00 3.259e+03 -0.977
## Categoryopioid_instrument_hand 1.095e+00 1.464e+00 3.259e+03 0.748
## Time:VisitV2 -2.181e-02 8.055e-03 3.259e+03 -2.707
## Age:SexMale -2.342e-01 6.045e-01 2.300e+01 -0.387
## InjectionUserTRUE:Categoryopioid_face_activities 7.640e-01 1.695e+00 3.259e+03 0.451
## InjectionUserTRUE:Categoryopioid_hand 8.393e-01 1.958e+00 3.259e+03 0.429
## InjectionUserTRUE:Categoryopioid_injection_hand 1.283e+00 1.958e+00 3.259e+03 0.655
## InjectionUserTRUE:Categoryopioid_injection_instrument 3.275e+00 1.495e+00 3.259e+03 2.190
## InjectionUserTRUE:Categoryopioid_instrument_hand 1.285e+00 1.696e+00 3.259e+03 0.758
## Pr(>|t|)
## (Intercept) 9.31e-06 ***
## Time 4.48e-06 ***
## VisitV2 0.00246 **
## Age 0.16583
## SexMale 0.84374
## InjectionUserTRUE 0.47502
## Categoryopioid_face_activities 0.26977
## Categoryopioid_hand 0.24393
## Categoryopioid_injection_hand 0.28015
## Categoryopioid_injection_instrument 0.32881
## Categoryopioid_instrument_hand 0.45476
## Time:VisitV2 0.00682 **
## Age:SexMale 0.70205
## InjectionUserTRUE:Categoryopioid_face_activities 0.65228
## InjectionUserTRUE:Categoryopioid_hand 0.66812
## InjectionUserTRUE:Categoryopioid_injection_hand 0.51222
## InjectionUserTRUE:Categoryopioid_injection_instrument 0.02857 *
## InjectionUserTRUE:Categoryopioid_instrument_hand 0.44862
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 12568.7 12568.7 1 3259.2 83.8112 < 2.2e-16 ***
## Visit 1377.2 1377.2 1 3259.9 9.1832 0.002462 **
## Age 285.9 285.9 1 23.0 1.9061 0.180666
## Sex 6.0 6.0 1 23.0 0.0397 0.843737
## InjectionUser 39.3 39.3 1 23.1 0.2617 0.613788
## Category 3445.4 689.1 5 3259.0 4.5949 0.000352 ***
## Time:Visit 1099.2 1099.2 1 3259.3 7.3300 0.006817 **
## Age:Sex 22.5 22.5 1 23.0 0.1501 0.702047
## InjectionUser:Category 917.3 183.5 5 3259.0 1.2233 0.295304
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('opioid'),], 'craving')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 26704
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.5357 -0.1324 0.1776 0.4725 3.4332
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 133.1 11.54
## Residual 186.3 13.65
## Number of obs: 3300, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 7.700e+01 1.387e+01 2.336e+01 5.551
## Time -3.645e-02 6.300e-03 3.259e+03 -5.786
## VisitV2 -3.165e+00 9.537e-01 3.260e+03 -3.319
## Age 4.851e-01 3.566e-01 2.300e+01 1.360
## SexMale -2.330e+00 2.300e+01 2.300e+01 -0.101
## InjectionUserTRUE -1.829e+00 5.525e+00 2.542e+01 -0.331
## Categoryopioid_face_activities 2.885e+00 1.632e+00 3.259e+03 1.768
## Categoryopioid_hand -1.747e+00 1.884e+00 3.259e+03 -0.927
## Categoryopioid_injection_hand 1.484e+00 1.884e+00 3.259e+03 0.788
## Categoryopioid_injection_instrument -1.286e+00 1.439e+00 3.259e+03 -0.894
## Categoryopioid_instrument_hand 1.925e+00 1.632e+00 3.259e+03 1.180
## Time:VisitV2 -1.274e-02 8.978e-03 3.259e+03 -1.419
## Age:SexMale -3.938e-02 5.981e-01 2.300e+01 -0.066
## InjectionUserTRUE:Categoryopioid_face_activities 3.055e-01 1.890e+00 3.259e+03 0.162
## InjectionUserTRUE:Categoryopioid_hand 1.123e+00 2.182e+00 3.259e+03 0.515
## InjectionUserTRUE:Categoryopioid_injection_hand 3.472e+00 2.182e+00 3.259e+03 1.591
## InjectionUserTRUE:Categoryopioid_injection_instrument 3.516e+00 1.666e+00 3.259e+03 2.110
## InjectionUserTRUE:Categoryopioid_instrument_hand 1.891e+00 1.890e+00 3.259e+03 1.000
## Pr(>|t|)
## (Intercept) 1.14e-05 ***
## Time 7.87e-09 ***
## VisitV2 0.000914 ***
## Age 0.186878
## SexMale 0.920165
## InjectionUserTRUE 0.743330
## Categoryopioid_face_activities 0.077145 .
## Categoryopioid_hand 0.353876
## Categoryopioid_injection_hand 0.430899
## Categoryopioid_injection_instrument 0.371506
## Categoryopioid_instrument_hand 0.238175
## Time:VisitV2 0.156081
## Age:SexMale 0.948075
## InjectionUserTRUE:Categoryopioid_face_activities 0.871566
## InjectionUserTRUE:Categoryopioid_hand 0.606749
## InjectionUserTRUE:Categoryopioid_injection_hand 0.111678
## InjectionUserTRUE:Categoryopioid_injection_instrument 0.034932 *
## InjectionUserTRUE:Categoryopioid_instrument_hand 0.317242
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 16945.9 16945.9 1 3259.3 90.9614 < 2.2e-16 ***
## Visit 2051.9 2051.9 1 3260.1 11.0143 0.0009141 ***
## Age 494.3 494.3 1 23.0 2.6530 0.1169728
## Sex 1.9 1.9 1 23.0 0.0103 0.9201647
## InjectionUser 0.1 0.1 1 23.1 0.0004 0.9837537
## Category 5978.7 1195.7 5 3259.0 6.4185 6.078e-06 ***
## Time:Visit 375.0 375.0 1 3259.4 2.0127 0.1560811
## Age:Sex 0.8 0.8 1 23.0 0.0043 0.9480752
## InjectionUser:Category 1324.8 265.0 5 3259.0 1.4223 0.2128160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lm_imagechar <- function(dset, characteristic){
this_lm <- lm(as.formula(paste(characteristic, '~ ImageSet')),
data = dset)
print(plot_model(this_lm, title = characteristic, show.values = TRUE, type = 'est'))
print(summary(this_lm))
print(anova(this_lm))
}
summarize_one_lm_imagechar(merged_summary_values, 'hue_mean')##
## Call:
## lm(formula = as.formula(paste(characteristic, "~ ImageSet")),
## data = dset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27365 -0.14219 -0.04265 0.10681 0.59849
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.19051 0.01493 12.757 < 2e-16 ***
## ImageSetmeth 0.08313 0.02112 3.936 9.85e-05 ***
## ImageSetopioid 0.06968 0.02112 3.299 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1683 on 378 degrees of freedom
## Multiple R-squared: 0.04511, Adjusted R-squared: 0.04005
## F-statistic: 8.928 on 2 and 378 DF, p-value: 0.0001628
##
## Analysis of Variance Table
##
## Response: hue_mean
## Df Sum Sq Mean Sq F value Pr(>F)
## ImageSet 2 0.5058 0.252875 8.9277 0.0001628 ***
## Residuals 378 10.7068 0.028325
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lm_imagechar(merged_summary_values, 'saturation_mean')##
## Call:
## lm(formula = as.formula(paste(characteristic, "~ ImageSet")),
## data = dset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30227 -0.12424 -0.01327 0.10830 0.44876
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.18970 0.01446 13.117 < 2e-16 ***
## ImageSetmeth 0.11257 0.02045 5.504 6.86e-08 ***
## ImageSetopioid 0.10954 0.02045 5.356 1.48e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.163 on 378 degrees of freedom
## Multiple R-squared: 0.09425, Adjusted R-squared: 0.08946
## F-statistic: 19.67 on 2 and 378 DF, p-value: 7.493e-09
##
## Analysis of Variance Table
##
## Response: saturation_mean
## Df Sum Sq Mean Sq F value Pr(>F)
## ImageSet 2 1.0448 0.52240 19.667 7.493e-09 ***
## Residuals 378 10.0406 0.02656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lm_imagechar(merged_summary_values, 'value_mean')##
## Call:
## lm(formula = as.formula(paste(characteristic, "~ ImageSet")),
## data = dset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42727 -0.15408 0.00472 0.10972 0.59092
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.45227 0.01884 24.004 <2e-16 ***
## ImageSetmeth -0.05318 0.02665 -1.996 0.0467 *
## ImageSetopioid -0.02398 0.02665 -0.900 0.3687
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2123 on 378 degrees of freedom
## Multiple R-squared: 0.01046, Adjusted R-squared: 0.005226
## F-statistic: 1.998 on 2 and 378 DF, p-value: 0.137
##
## Analysis of Variance Table
##
## Response: value_mean
## Df Sum Sq Mean Sq F value Pr(>F)
## ImageSet 2 0.1802 0.090085 1.9982 0.137
## Residuals 378 17.0419 0.045084
#ddq_data <- read.csv('DDQ_Scores.csv')
ddq_data <- read.csv('DCR_DDQScores-anonymized.csv')
summarize_one_lme_ddq <- function(dset, subscore){
this_lme <- lmer(Value ~ Time + (1 | id),
data = dset[dset$Subscore == subscore,])
print(plot_model(this_lme, title = subscore, show.values = TRUE, type = 'std'))
print(summary(this_lme))
print(anova(this_lme))
}
summarize_one_lme_ddq(ddq_data, 'Desire')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time + (1 | id)
## Data: dset[dset$Subscore == subscore, ]
##
## REML criterion at convergence: 215.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.94412 -0.25263 -0.06047 0.21956 3.10401
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.6939 0.8330
## Residual 0.1921 0.4383
## Number of obs: 112, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.70918 0.17789 38.02770 9.608 1.02e-11 ***
## TimeT2 0.03571 0.11715 81.00000 0.305 0.761
## TimeT3 -0.07143 0.11715 81.00000 -0.610 0.544
## TimeT4 -0.08673 0.11715 81.00000 -0.740 0.461
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TimeT2 TimeT3
## TimeT2 -0.329
## TimeT3 -0.329 0.500
## TimeT4 -0.329 0.500 0.500
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 0.28426 0.094752 3 81 0.4932 0.688
summarize_one_lme_ddq(ddq_data, 'Neg')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time + (1 | id)
## Data: dset[dset$Subscore == subscore, ]
##
## REML criterion at convergence: 338.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3745 -0.3277 -0.0669 0.2760 2.3180
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 2.5480 1.5963
## Residual 0.5674 0.7533
## Number of obs: 112, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.3929 0.3336 35.9197 7.174 1.97e-08 ***
## TimeT2 -0.1339 0.2013 81.0000 -0.665 0.5078
## TimeT3 -0.1964 0.2013 81.0000 -0.976 0.3321
## TimeT4 -0.4018 0.2013 81.0000 -1.996 0.0493 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TimeT2 TimeT3
## TimeT2 -0.302
## TimeT3 -0.302 0.500
## TimeT4 -0.302 0.500 0.500
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 2.3504 0.78348 3 81 1.3807 0.2546
summarize_one_lme_ddq(ddq_data, 'Control')## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time + (1 | id)
## Data: dset[dset$Subscore == subscore, ]
##
## REML criterion at convergence: 332.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06379 -0.36536 -0.07568 0.27796 2.94957
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 2.3385 1.5292
## Residual 0.5472 0.7397
## Number of obs: 112, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.14286 0.32103 36.36237 6.675 8.38e-08 ***
## TimeT2 0.16071 0.19770 81.00000 0.813 0.419
## TimeT3 -0.05357 0.19770 81.00000 -0.271 0.787
## TimeT4 -0.30357 0.19770 81.00000 -1.535 0.129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TimeT2 TimeT3
## TimeT2 -0.308
## TimeT3 -0.308 0.500
## TimeT4 -0.308 0.500 0.500
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 3.1138 1.0379 3 81 1.8968 0.1367
need: all images craving|arousal|valence|absvalence ~ drug (M/O/N) + age + gender
library(multcomp)## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:sm':
##
## muscle
## The following object is masked from 'package:dplyr':
##
## select
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
## The following object is masked from 'package:sm':
##
## geyser
summarize_one_lme_allimages <- function(dset, RatingType){
this_lme <- lmer(Value ~ ImageSet + Age + Sex + (1 | id),
data = dset[dset$RatingType == RatingType,])
print(paste("LME OF ALL IMAGES FOR:", RatingType))
print(plot_model(this_lme, title = paste(RatingType), show.values = TRUE, type = 'std', order.terms = c(3,4,1,2)))
print(summary(this_lme))
print(anova(this_lme))
this_lme.ht <- glht(this_lme, linfct = mcp(ImageSet = "Tukey"))
summary(this_lme.ht, test=adjusted("none"))
summary(this_lme.ht)
}
summarize_one_lme_allimages(subject_data, 'craving')## [1] "LME OF ALL IMAGES FOR: craving"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ ImageSet + Age + Sex + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 86227.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0716 -0.6391 0.2088 0.4583 5.3166
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 32.66 5.715
## Residual 351.75 18.755
## Number of obs: 9900, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.9670 5.4414 25.1223 0.729 0.473
## ImageSetmeth 76.1436 0.4617 9870.0491 164.914 <2e-16 ***
## ImageSetopioid 76.4842 0.4617 9870.0491 165.652 <2e-16 ***
## Age 0.2087 0.1376 25.0025 1.517 0.142
## SexMale -0.4920 2.2159 25.0336 -0.222 0.826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ImgStm ImgStp Age
## ImageSetmth -0.042
## ImageSetopd -0.042 0.500
## Age -0.950 0.000 0.000
## SexMale -0.219 0.000 0.000 -0.015
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ImageSet 12812590 6406295 2 9870 18212.5694 <2e-16 ***
## Age 809 809 1 25 2.3011 0.1418
## Sex 17 17 1 25 0.0493 0.8261
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = Value ~ ImageSet + Age + Sex + (1 | id), data = dset[dset$RatingType ==
## RatingType, ])
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## meth - control == 0 76.1436 0.4617 164.914 <1e-04 ***
## opioid - control == 0 76.4842 0.4617 165.652 <1e-04 ***
## opioid - meth == 0 0.3406 0.4617 0.738 0.741
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summarize_one_lme_allimages(subject_data, 'arousal')## [1] "LME OF ALL IMAGES FOR: arousal"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ ImageSet + Age + Sex + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 36788.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9768 -0.5700 0.0353 0.7985 4.4307
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 2.793 1.671
## Residual 2.364 1.537
## Number of obs: 9900, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.340e-01 1.568e+00 2.501e+01 0.404 0.6894
## ImageSetmeth 2.169e+00 3.785e-02 9.870e+03 57.310 <2e-16 ***
## ImageSetopioid 2.165e+00 3.785e-02 9.870e+03 57.198 <2e-16 ***
## Age 1.560e-02 3.968e-02 2.500e+01 0.393 0.6975
## SexMale 2.218e+00 6.390e-01 2.500e+01 3.471 0.0019 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ImgStm ImgStp Age
## ImageSetmth -0.012
## ImageSetopd -0.012 0.500
## Age -0.951 0.000 0.000
## SexMale -0.219 0.000 0.000 -0.014
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ImageSet 10330.7 5165.3 2 9870 2185.3590 < 2.2e-16 ***
## Age 0.4 0.4 1 25 0.1546 0.697523
## Sex 28.5 28.5 1 25 12.0468 0.001899 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = Value ~ ImageSet + Age + Sex + (1 | id), data = dset[dset$RatingType ==
## RatingType, ])
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## meth - control == 0 2.169091 0.037848 57.310 <1e-06 ***
## opioid - control == 0 2.164848 0.037848 57.198 <1e-06 ***
## opioid - meth == 0 -0.004242 0.037848 -0.112 0.993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summarize_one_lme_allimages(subject_data, 'valence')## [1] "LME OF ALL IMAGES FOR: valence"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ ImageSet + Age + Sex + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 34641.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2479 -0.4175 0.0989 0.4949 4.1855
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.038 1.019
## Residual 1.908 1.381
## Number of obs: 9897, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.12136 0.95728 25.02159 3.261 0.0032 **
## ImageSetmeth -0.44533 0.03402 9867.00911 -13.092 <2e-16 ***
## ImageSetopioid -0.44667 0.03401 9867.00900 -13.134 <2e-16 ***
## Age 0.03840 0.02423 25.00064 1.585 0.1256
## SexMale 0.53025 0.39011 25.00638 1.359 0.1862
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ImgStm ImgStp Age
## ImageSetmth -0.018
## ImageSetopd -0.018 0.500
## Age -0.951 0.000 0.000
## SexMale -0.219 0.000 0.000 -0.014
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ImageSet 437.55 218.773 2 9867 114.6380 <2e-16 ***
## Age 4.79 4.794 1 25 2.5120 0.1256
## Sex 3.53 3.526 1 25 1.8475 0.1862
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = Value ~ ImageSet + Age + Sex + (1 | id), data = dset[dset$RatingType ==
## RatingType, ])
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## meth - control == 0 -0.445328 0.034017 -13.092 <1e-07 ***
## opioid - control == 0 -0.446667 0.034009 -13.134 <1e-07 ***
## opioid - meth == 0 -0.001339 0.034017 -0.039 0.999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summarize_one_lme_allimages(subject_data, 'valence_fromneutral')## [1] "LME OF ALL IMAGES FOR: valence_fromneutral"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ ImageSet + Age + Sex + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 30388.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9879 -0.7091 -0.2381 0.4710 4.0593
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.7221 0.8498
## Residual 1.2413 1.1141
## Number of obs: 9897, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.39794 0.79833 25.01834 -0.498 0.623
## ImageSetmeth 0.83280 0.02743 9867.00669 30.357 <2e-16 ***
## ImageSetopioid 0.78788 0.02743 9867.00658 28.726 <2e-16 ***
## Age 0.02146 0.02020 24.99875 1.062 0.298
## SexMale 0.20780 0.32534 25.00413 0.639 0.529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ImgStm ImgStp Age
## ImageSetmth -0.017
## ImageSetopd -0.017 0.500
## Age -0.951 0.000 0.000
## SexMale -0.219 0.000 0.000 -0.014
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ImageSet 1447.71 723.86 2 9867 583.1594 <2e-16 ***
## Age 1.40 1.40 1 25 1.1277 0.2984
## Sex 0.51 0.51 1 25 0.4080 0.5288
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = Value ~ ImageSet + Age + Sex + (1 | id), data = dset[dset$RatingType ==
## RatingType, ])
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## meth - control == 0 0.83280 0.02743 30.357 <1e-04 ***
## opioid - control == 0 0.78788 0.02743 28.726 <1e-04 ***
## opioid - meth == 0 -0.04493 0.02743 -1.638 0.23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
drug images: craving|methToOpioid ~ cat drug + time * session
Need category to have the same values for opioid/meth here–e.g. instrument_hand instead of opioid_instrument_hand
#meth_instrument opioid_injection_instrument
#meth_hand opioid_hand
#meth_injection_hand opioid_injection_hand
#meth_instrument_hand opioid_instrument_hand
#meth_face_activities opioid_face_activities
#meth opioid
subject_data$UnifiedCategory <- NA
subject_data$UnifiedCategory[subject_data$Category %in% c('meth', 'opioid')] <- 'drug'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_instrument', 'opioid_injection_instrument')] <- 'instrument'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_hand', 'opioid_hand')] <- 'drug_hand'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_injection_hand', 'opioid_injection_hand')] <- 'injection_hand'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_instrument_hand', 'opioid_instrument_hand')] <- 'instrument_hand'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_face_activities', 'opioid_face_activities')] <- 'face_activities'
subject_data$MethToOpioid <- NA
subject_data$MethToOpioid[subject_data$RelatedCategory == 'Meth'] <- -1
subject_data$MethToOpioid[subject_data$RelatedCategory == 'Neither'] <- 0
subject_data$MethToOpioid[subject_data$RelatedCategory == 'Opioid'] <- 1
subject_data$MethToOpioid[subject_data$RelatedCategory == 'Both'] <- 0
summarize_one_lme_drugcats <- function(dset, RatingType){
this_lme <- lmer(Value ~ UnifiedCategory * ImageSet + Time * Visit + (1 | id),
data = dset[dset$RatingType == RatingType,])
print(paste("LME OF Drug Images FOR:", RatingType))
print(plot_model(this_lme, title = paste(RatingType), show.values = TRUE, type = 'std'))
print(summary(this_lme))
print(anova(this_lme))
}
summarize_one_lme_drugcats(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'craving')## [1] "LME OF Drug Images FOR: craving"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ UnifiedCategory * ImageSet + Time * Visit + (1 | id)
## Data: dset[dset$RatingType == RatingType, ]
##
## REML criterion at convergence: 53644.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.4418 -0.1264 0.1650 0.4942 3.4449
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 132.5 11.51
## Residual 194.5 13.95
## Number of obs: 6600, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.420e+01 2.299e+00 3.329e+01 40.973 < 2e-16 ***
## UnifiedCategorydrug_hand -1.988e+00 9.711e-01 6.558e+03 -2.048 0.040634 *
## UnifiedCategoryface_activities -1.302e+00 8.412e-01 6.558e+03 -1.548 0.121746
## UnifiedCategoryinjection_hand -6.386e-01 9.714e-01 6.558e+03 -0.657 0.510990
## UnifiedCategoryinstrument -2.950e+00 7.422e-01 6.558e+03 -3.975 7.13e-05 ***
## UnifiedCategoryinstrument_hand -2.969e+00 8.415e-01 6.558e+03 -3.528 0.000422 ***
## ImageSetopioid -3.334e+00 8.411e-01 6.558e+03 -3.964 7.45e-05 ***
## Time -2.831e-02 4.519e-03 6.558e+03 -6.264 4.00e-10 ***
## VisitV2 -1.484e+00 6.910e-01 6.559e+03 -2.148 0.031770 *
## UnifiedCategorydrug_hand:ImageSetopioid 1.079e+00 1.373e+00 6.558e+03 0.786 0.431959
## UnifiedCategoryface_activities:ImageSetopioid 4.394e+00 1.189e+00 6.558e+03 3.694 0.000222 ***
## UnifiedCategoryinjection_hand:ImageSetopioid 4.739e+00 1.374e+00 6.558e+03 3.450 0.000564 ***
## UnifiedCategoryinstrument:ImageSetopioid 4.280e+00 1.049e+00 6.558e+03 4.081 4.55e-05 ***
## UnifiedCategoryinstrument_hand:ImageSetopioid 6.292e+00 1.190e+00 6.558e+03 5.289 1.27e-07 ***
## Time:VisitV2 -2.774e-02 6.474e-03 6.558e+03 -4.284 1.86e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## UnifiedCategory 5567 1113 5 6558.0 5.7243 2.821e-05 ***
## ImageSet 24 24 1 6558.0 0.1245 0.72422
## Time 32938 32938 1 6558.2 169.3537 < 2.2e-16 ***
## Visit 897 897 1 6559.0 4.6128 0.03177 *
## UnifiedCategory:ImageSet 7214 1443 5 6558.0 7.4185 5.999e-07 ***
## Time:Visit 3569 3569 1 6558.1 18.3516 1.863e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dset <- subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),]
this_lme <- lmer(MethToOpioid ~ UnifiedCategory * ImageSet + Time * Visit + (1 | id),
data = dset)
print(paste("LME OF Drug Images FOR:", 'MethToOpioid'))## [1] "LME OF Drug Images FOR: MethToOpioid"
print(plot_model(this_lme, title = paste('MethToOpioid'), show.values = TRUE, type = 'std'))print(summary(this_lme))## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MethToOpioid ~ UnifiedCategory * ImageSet + Time * Visit + (1 | id)
## Data: dset
##
## REML criterion at convergence: 11106.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3056 -0.6445 -0.2287 0.8722 3.7004
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.02873 0.1695
## Residual 0.30705 0.5541
## Number of obs: 6597, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.033e+00 4.361e-02 8.478e+01 -23.693 < 2e-16 ***
## UnifiedCategorydrug_hand 3.599e-01 3.862e-02 6.555e+03 9.319 < 2e-16 ***
## UnifiedCategoryface_activities 3.032e-01 3.342e-02 6.555e+03 9.072 < 2e-16 ***
## UnifiedCategoryinjection_hand 8.629e-01 3.864e-02 6.555e+03 22.336 < 2e-16 ***
## UnifiedCategoryinstrument 3.153e-01 2.949e-02 6.555e+03 10.690 < 2e-16 ***
## UnifiedCategoryinstrument_hand 1.660e-01 3.343e-02 6.555e+03 4.965 7.05e-07 ***
## ImageSetopioid 1.585e+00 3.342e-02 6.555e+03 47.415 < 2e-16 ***
## Time 4.816e-04 1.796e-04 6.556e+03 2.682 0.00733 **
## VisitV2 1.357e-01 2.745e-02 6.561e+03 4.943 7.89e-07 ***
## UnifiedCategorydrug_hand:ImageSetopioid -5.983e-01 5.459e-02 6.555e+03 -10.959 < 2e-16 ***
## UnifiedCategoryface_activities:ImageSetopioid -7.466e-01 4.726e-02 6.555e+03 -15.796 < 2e-16 ***
## UnifiedCategoryinjection_hand:ImageSetopioid -1.336e+00 5.460e-02 6.555e+03 -24.466 < 2e-16 ***
## UnifiedCategoryinstrument:ImageSetopioid -6.712e-01 4.168e-02 6.555e+03 -16.101 < 2e-16 ***
## UnifiedCategoryinstrument_hand:ImageSetopioid -5.402e-01 4.726e-02 6.555e+03 -11.430 < 2e-16 ***
## Time:VisitV2 -3.814e-04 2.573e-04 6.556e+03 -1.483 0.13821
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(summary(this_lme), correlation=TRUE) or
## vcov(summary(this_lme)) if you need it
print(anova(this_lme))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## UnifiedCategory 45.01 9.00 5 6555.1 29.3151 < 2.2e-16 ***
## ImageSet 1257.39 1257.39 1 6555.1 4095.0401 < 2.2e-16 ***
## Time 1.57 1.57 1 6556.1 5.1050 0.02389 *
## Visit 7.50 7.50 1 6561.4 24.4316 7.891e-07 ***
## UnifiedCategory:ImageSet 197.43 39.49 5 6555.1 128.5940 < 2.2e-16 ***
## Time:Visit 0.67 0.67 1 6556.0 2.1983 0.13821
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_time <- function(dset, formula_to_use){
this_lme <- lmer(formula_to_use,
data = dset)
print(paste("LME OF ALL IMAGES FOR: Time On Page"))
print(plot_model(this_lme, title = paste('Time On Page for Craving Ratings'), show.values = TRUE, type = 'std', order.terms = c(3,4,1,2)))
print(summary(this_lme))
print(anova(this_lme))
this_lme.ht <- glht(this_lme, linfct = mcp(ImageSet = "Tukey"))
summary(this_lme.ht, test=adjusted("none"))
summary(this_lme.ht)
}
craving_data <- subject_data[subject_data$RatingType == 'craving',]
data_touse <- craving_data[craving_data$millisecondsOnPage < 20000,]
print(paste0("Fraction Excluded for >20s on page: ", 1 - (nrow(data_touse) / nrow(craving_data))))## [1] "Fraction Excluded for >20s on page: 0.0060606060606061"
summarize_one_lme_time(data_touse, as.formula('millisecondsOnPage / 1000 ~ ImageSet + Age + Sex +(1 | id)'))## [1] "LME OF ALL IMAGES FOR: Time On Page"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
## Data: dset
##
## REML criterion at convergence: 39808.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1038 -0.5020 -0.2116 0.1848 9.5495
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.7481 0.865
## Residual 3.3002 1.817
## Number of obs: 9840, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.462e+00 8.159e-01 2.504e+01 1.791 0.085348 .
## ImageSetmeth 1.745e-01 4.485e-02 9.810e+03 3.890 0.000101 ***
## ImageSetopioid 3.084e-01 4.486e-02 9.810e+03 6.874 6.62e-12 ***
## Age 3.856e-02 2.064e-02 2.499e+01 1.868 0.073535 .
## SexMale 3.581e-01 3.324e-01 2.501e+01 1.077 0.291712
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ImgStm ImgStp Age
## ImageSetmth -0.027
## ImageSetopd -0.027 0.500
## Age -0.951 0.000 0.000
## SexMale -0.219 0.000 0.000 -0.014
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ImageSet 156.853 78.426 2 9810 23.7642 5.061e-11 ***
## Age 11.515 11.515 1 25 3.4893 0.07353 .
## Sex 3.829 3.829 1 25 1.1602 0.29171
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = formula_to_use, data = dset)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## meth - control == 0 0.17449 0.04485 3.890 < 0.001 ***
## opioid - control == 0 0.30836 0.04486 6.874 < 0.001 ***
## opioid - meth == 0 0.13387 0.04487 2.984 0.00803 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summarize_one_lme_time(data_touse, as.formula('millisecondsOnPage / 1000 ~ ImageSet + Age + Sex + Time * Visit + (1 | id)'))## [1] "LME OF ALL IMAGES FOR: Time On Page"
## Number of values in `order.terms` does not match number of terms. Terms are not sorted.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
## Data: dset
##
## REML criterion at convergence: 38644
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2602 -0.5112 -0.1761 0.2246 10.2432
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.7705 0.8778
## Residual 2.9220 1.7094
## Number of obs: 9840, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.825e+00 8.284e-01 2.517e+01 3.410 0.0022 **
## ImageSetmeth 2.529e-01 4.231e-02 9.807e+03 5.979 2.33e-09 ***
## ImageSetopioid 3.763e-01 4.229e-02 9.807e+03 8.897 < 2e-16 ***
## Age 3.909e-02 2.093e-02 2.498e+01 1.868 0.0736 .
## SexMale 3.364e-01 3.371e-01 2.499e+01 0.998 0.3279
## Time -1.180e-02 4.664e-04 9.807e+03 -25.305 < 2e-16 ***
## VisitV2 -1.240e+00 6.900e-02 9.809e+03 -17.967 < 2e-16 ***
## Time:VisitV2 5.510e-03 6.642e-04 9.807e+03 8.296 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ImgStm ImgStp Age SexMal Time VistV2
## ImageSetmth -0.023
## ImageSetopd -0.024 0.502
## Age -0.950 0.000 0.000
## SexMale -0.219 0.000 0.000 -0.014
## Time -0.049 -0.046 -0.036 0.000 0.000
## VisitV2 -0.041 0.005 0.010 -0.001 0.001 0.605
## Time:VistV2 0.036 -0.005 -0.012 0.000 0.000 -0.700 -0.864
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ImageSet 240.22 120.11 2 9807.0 41.105 < 2e-16 ***
## Age 10.19 10.19 1 25.0 3.488 0.07359 .
## Sex 2.91 2.91 1 25.0 0.996 0.32785
## Time 2155.99 2155.99 1 9807.0 737.839 < 2e-16 ***
## Visit 943.31 943.31 1 9808.8 322.827 < 2e-16 ***
## Time:Visit 201.13 201.13 1 9807.0 68.831 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = formula_to_use, data = dset)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## meth - control == 0 0.25294 0.04231 5.979 < 0.001 ***
## opioid - control == 0 0.37628 0.04229 8.897 < 0.001 ***
## opioid - meth == 0 0.12333 0.04222 2.921 0.00976 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
total_values <- aggregate(subject_data$millisecondsOnPage,
by = list(id = subject_data$id, Visit = subject_data$Visit),
FUN = sum, na.rm = TRUE)
names(total_values) <- c('id', 'Visit', 'TaskCompletionTime')
#convert to seconds
total_values$TaskCompletionTime <- total_values$TaskCompletionTime / 1000
print("Total task completion time for each visit:")## [1] "Total task completion time for each visit:"
CreateTableOne(vars = 'TaskCompletionTime', data = total_values, strata = 'Visit')## Stratified by Visit
## V1 V2 p test
## n 28 27
## TaskCompletionTime (mean (sd)) 4793.03 (1653.44) 3583.79 (1081.61) 0.002
image_total_values <- aggregate(subject_data$millisecondsOnPage,
by = list(id = subject_data$id, ImageSetFile = subject_data$ImageSetFile),
FUN = sum, na.rm = TRUE)
names(image_total_values) <- c('id', 'ImageSetFile', 'ImageCompletionTime')
#convert to seconds
image_total_values$ImageCompletionTime <- image_total_values$ImageCompletionTime / 1000
subject_data <- merge(subject_data, image_total_values)
summarize_one_lme_imagetime <- function(dset, formula_to_use){
this_lme <- lmer(formula_to_use,
data = dset)
print(paste("LME OF DRUG IMAGES FOR: Time On Image"))
print(plot_model(this_lme, title = paste('Time On Image for Craving Ratings'), show.values = TRUE, type = 'std', order.terms = c(3,4,2,1)))
print(summary(this_lme))
print(anova(this_lme))
#this_lme.ht <- glht(this_lme, linfct = mcp(ImageSet = "Tukey"))
#summary(this_lme.ht, test=adjusted("none"))
#summary(this_lme.ht)
}
craving_data <- subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet %in% c('meth', 'opioid'),]
names(craving_data)[names(craving_data) == 'Value'] <- 'craving'
data_touse <- craving_data[craving_data$ImageCompletionTime < 100,]
print(paste0("Fraction Excluded for >100s on image: ", 1 - (nrow(data_touse) / nrow(craving_data))))## [1] "Fraction Excluded for >100s on image: 0.0106060606060606"
summarize_one_lme_imagetime(data_touse, as.formula('ImageCompletionTime ~ craving + ImageSet + Time * Visit + (1 | id)'))## [1] "LME OF DRUG IMAGES FOR: Time On Image"
## Number of values in `order.terms` does not match number of terms. Terms are not sorted.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
## Data: dset
##
## REML criterion at convergence: 48833.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2938 -0.5597 -0.1788 0.2675 8.5762
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 36.87 6.072
## Residual 101.32 10.066
## Number of obs: 6530, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.965e+01 1.459e+00 6.835e+01 20.315 < 2e-16 ***
## craving 3.829e-02 8.876e-03 6.506e+03 4.314 1.63e-05 ***
## ImageSetopioid 1.134e+00 2.492e-01 6.497e+03 4.552 5.40e-06 ***
## Time -9.147e-02 3.308e-03 6.498e+03 -27.651 < 2e-16 ***
## VisitV2 -7.773e+00 5.042e-01 6.499e+03 -15.415 < 2e-16 ***
## Time:VisitV2 3.021e-02 4.715e-03 6.498e+03 6.408 1.58e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cravng ImgStp Time VistV2
## craving -0.562
## ImageSetopd -0.082 -0.010
## Time -0.252 0.077 0.010
## VisitV2 -0.184 0.026 0.007 0.607
## Time:VistV2 0.117 0.053 -0.009 -0.693 -0.864
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## craving 1885 1885 1 6506.3 18.608 1.630e-05 ***
## ImageSet 2100 2100 1 6497.0 20.725 5.400e-06 ***
## Time 103901 103901 1 6499.5 1025.490 < 2.2e-16 ***
## Visit 24075 24075 1 6499.0 237.618 < 2.2e-16 ***
## Time:Visit 4160 4160 1 6497.6 41.061 1.579e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_imagetime(data_touse, as.formula('ImageCompletionTime ~ craving + Time * Visit + (1 | id)'))## [1] "LME OF DRUG IMAGES FOR: Time On Image"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
## Data: dset
##
## REML criterion at convergence: 48853.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3473 -0.5595 -0.1793 0.2651 8.5086
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 36.86 6.071
## Residual 101.63 10.081
## Number of obs: 6530, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.019e+01 1.455e+00 6.761e+01 20.746 < 2e-16 ***
## craving 3.870e-02 8.889e-03 6.507e+03 4.353 1.36e-05 ***
## Time -9.162e-02 3.313e-03 6.499e+03 -27.656 < 2e-16 ***
## VisitV2 -7.789e+00 5.050e-01 6.500e+03 -15.423 < 2e-16 ***
## Time:VisitV2 3.041e-02 4.722e-03 6.499e+03 6.439 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cravng Time VistV2
## craving -0.566
## Time -0.253 0.077
## VisitV2 -0.184 0.026 0.607
## Time:VistV2 0.117 0.053 -0.693 -0.864
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## craving 1926 1926 1 6507.0 18.951 1.362e-05 ***
## Time 104049 104049 1 6500.5 1023.842 < 2.2e-16 ***
## Visit 24173 24173 1 6500.1 237.863 < 2.2e-16 ***
## Time:Visit 4213 4213 1 6498.6 41.459 1.290e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_imagetime(data_touse, as.formula('ImageCompletionTime ~ craving + (1 | id)'))## [1] "LME OF DRUG IMAGES FOR: Time On Image"
## Number of values in `order.terms` does not match number of terms. Terms are not sorted.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
## Data: dset
##
## REML criterion at convergence: 50129.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0339 -0.5448 -0.2364 0.2137 7.6031
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 35.35 5.946
## Residual 124.01 11.136
## Number of obs: 6530, groups: id, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.353e+01 1.411e+00 6.443e+01 9.589 5.03e-14 ***
## craving 1.043e-01 9.574e-03 6.489e+03 10.897 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## craving -0.597
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## craving 14724 14724 1 6489.4 118.73 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
this_lme <- lmer(millisecondsOnPage / 1000 ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit +(1 | id),
data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'meth',])## Warning: Some predictor variables are on very different scales: consider rescaling
## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'meth RT to Craving Ratings', show.values = TRUE, type = 'std', order.terms = c(1,2,8,9,10,3,4,6,5,7)))print(summary(this_lme))## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: millisecondsOnPage/1000 ~ Age + Sex + MethOnsetAge + MethMonthsUsed +
## MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit + (1 | id)
## Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==
## "meth", ]
##
## REML criterion at convergence: 17740.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.012 -0.293 -0.116 0.071 32.336
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.7159 0.8461
## Residual 14.9185 3.8625
## Number of obs: 3180, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.341e+00 1.529e+00 1.980e+01 2.840 0.0102 *
## Age 5.741e-02 2.813e-02 1.918e+01 2.040 0.0553 .
## SexMale -3.078e-01 4.682e-01 1.916e+01 -0.657 0.5188
## MethOnsetAge 2.185e-02 4.571e-02 1.922e+01 0.478 0.6381
## MethMonthsUsed -6.900e-04 2.857e-03 1.905e+01 -0.241 0.8118
## MethDollars 2.532e-06 9.471e-05 1.946e+01 0.027 0.9789
## MethDaysLastMonth -3.776e-02 1.845e-02 1.925e+01 -2.047 0.0546 .
## MethDaysAbstinant -1.063e-03 6.759e-04 1.908e+01 -1.573 0.1323
## Time -1.366e-02 1.804e-03 3.153e+03 -7.569 4.89e-14 ***
## VisitV2 -1.305e+00 2.778e-01 3.158e+03 -4.697 2.76e-06 ***
## Time:VisitV2 4.547e-03 2.591e-03 3.153e+03 1.755 0.0793 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Age SexMal MthOnA MthMnU MthDll MthDLM MthDyA Time VistV2
## Age -0.526
## SexMale -0.053 -0.007
## MethOnsetAg -0.570 -0.147 -0.390
## MthMnthsUsd -0.216 -0.386 -0.325 0.573
## MethDollars -0.341 0.121 -0.017 0.174 -0.058
## MthDysLstMn -0.145 -0.115 0.350 -0.023 -0.069 -0.438
## MthDysAbstn -0.504 0.076 0.329 0.055 0.098 0.390 0.053
## Time -0.107 -0.001 -0.007 0.006 0.007 -0.006 -0.003 -0.013
## VisitV2 -0.075 -0.012 0.006 -0.007 0.010 -0.010 0.003 -0.010 0.602
## Time:VistV2 0.071 0.006 0.000 -0.001 -0.006 -0.001 0.005 0.005 -0.695 -0.868
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Age 62.12 62.12 1 19.18 4.1636 0.05531 .
## Sex 6.45 6.45 1 19.16 0.4321 0.51880
## MethOnsetAge 3.41 3.41 1 19.22 0.2284 0.63808
## MethMonthsUsed 0.87 0.87 1 19.05 0.0583 0.81175
## MethDollars 0.01 0.01 1 19.46 0.0007 0.97894
## MethDaysLastMonth 62.50 62.50 1 19.25 4.1893 0.05458 .
## MethDaysAbstinant 36.89 36.89 1 19.08 2.4730 0.13225
## Time 1148.29 1148.29 1 3154.39 76.9707 < 2.2e-16 ***
## Visit 329.06 329.06 1 3158.00 22.0574 2.759e-06 ***
## Time:Visit 45.96 45.96 1 3152.70 3.0805 0.07934 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
this_lme <- lmer(millisecondsOnPage / 1000 ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed + OpioidOrHeroinDollars +
OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant + Time * Visit + (1 | id),
data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'opioid',])## Warning: Some predictor variables are on very different scales: consider rescaling
## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'opioid RT to Craving Ratings', show.values = TRUE, type = 'std', order.terms = c(1,2,8, 9, 10, 3,4,6,5,7)))print(summary(this_lme))## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: millisecondsOnPage/1000 ~ Age + Sex + OpioidOrHeroinOnsetAge +
## OpioidOrHeroinMonthsUsed + OpioidOrHeroinDollars + OpioidOrHeroinDaysLastMonth +
## OpioidOrHeroinDaysAbstinant + Time * Visit + (1 | id)
## Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==
## "opioid", ]
##
## REML criterion at convergence: 17831.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.026 -0.323 -0.127 0.084 36.956
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.056 1.028
## Residual 15.314 3.913
## Number of obs: 3180, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.250e+00 1.500e+00 1.939e+01 2.833 0.0105 *
## Age 4.577e-02 4.288e-02 1.891e+01 1.067 0.2993
## SexMale 2.286e-01 4.654e-01 1.912e+01 0.491 0.6289
## OpioidOrHeroinOnsetAge 1.486e-02 4.888e-02 1.920e+01 0.304 0.7644
## OpioidOrHeroinMonthsUsed -6.160e-04 3.490e-03 1.903e+01 -0.176 0.8618
## OpioidOrHeroinDollars -1.590e-04 1.643e-04 1.918e+01 -0.968 0.3451
## OpioidOrHeroinDaysLastMonth -1.926e-02 2.604e-02 1.887e+01 -0.740 0.4686
## OpioidOrHeroinDaysAbstinant -1.380e-04 1.879e-04 1.894e+01 -0.735 0.4715
## Time -1.509e-02 1.833e-03 3.153e+03 -8.235 2.61e-16 ***
## VisitV2 -1.478e+00 2.772e-01 3.157e+03 -5.332 1.04e-07 ***
## Time:VisitV2 5.449e-03 2.611e-03 3.152e+03 2.087 0.0370 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Age SexMal OpOHOA OpOHMU OpdOHD OOHDLM OpOHDA Time VistV2
## Age -0.750
## SexMale -0.065 0.021
## OpdOrHrnOnA 0.083 -0.616 -0.149
## OpdOrHrnMnU 0.125 -0.397 -0.336 0.409
## OpdOrHrnDll 0.058 0.119 0.215 -0.421 -0.250
## OpdOrHrnDLM -0.629 0.486 0.056 -0.261 -0.392 -0.157
## OpdOrHrnDyA -0.177 -0.045 -0.160 -0.003 0.227 0.016 0.236
## Time -0.104 -0.002 0.003 -0.006 -0.003 0.004 -0.005 0.001
## VisitV2 -0.080 -0.004 0.006 -0.007 -0.008 0.010 -0.005 -0.003 0.603
## Time:VistV2 0.069 0.008 0.002 -0.002 0.002 -0.002 0.005 -0.002 -0.702 -0.863
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Age 17.44 17.44 1 18.91 1.1390 0.29930
## Sex 3.70 3.70 1 19.12 0.2413 0.62886
## OpioidOrHeroinOnsetAge 1.41 1.41 1 19.20 0.0924 0.76444
## OpioidOrHeroinMonthsUsed 0.48 0.48 1 19.03 0.0311 0.86177
## OpioidOrHeroinDollars 14.35 14.35 1 19.18 0.9371 0.34509
## OpioidOrHeroinDaysLastMonth 8.38 8.38 1 18.87 0.5471 0.46861
## OpioidOrHeroinDaysAbstinant 8.27 8.27 1 18.94 0.5399 0.47146
## Time 1374.11 1374.11 1 3152.11 89.7303 < 2.2e-16 ***
## Visit 435.33 435.33 1 3156.68 28.4271 1.041e-07 ***
## Time:Visit 66.69 66.69 1 3152.34 4.3547 0.03699 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
this_lme <- lmer( ImageCompletionTime ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit +(1 | id),
data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'meth',])## Warning: Some predictor variables are on very different scales: consider rescaling
## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'meth Image Completion Time', show.values = TRUE, type = 'std', order.terms = c(1,2,8,9,10,3,4,6,5,7)))print(summary(this_lme))## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: ImageCompletionTime ~ Age + Sex + MethOnsetAge + MethMonthsUsed +
## MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit + (1 | id)
## Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==
## "meth", ]
##
## REML criterion at convergence: 30424.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8646 -0.2624 -0.1015 0.0871 31.3912
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 48.59 6.971
## Residual 815.56 28.558
## Number of obs: 3180, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.277e+01 1.239e+01 1.957e+01 1.838 0.0813 .
## Age 3.293e-01 2.283e-01 1.906e+01 1.442 0.1655
## SexMale -2.542e+00 3.800e+00 1.904e+01 -0.669 0.5116
## MethOnsetAge 3.369e-01 3.709e-01 1.909e+01 0.908 0.3751
## MethMonthsUsed 7.553e-03 2.320e-02 1.895e+01 0.326 0.7483
## MethDollars -8.147e-05 7.681e-04 1.930e+01 -0.106 0.9166
## MethDaysLastMonth -3.635e-02 1.497e-01 1.911e+01 -0.243 0.8108
## MethDaysAbstinant -3.952e-03 5.487e-03 1.897e+01 -0.720 0.4801
## Time -1.059e-01 1.334e-02 3.153e+03 -7.943 2.73e-15 ***
## VisitV2 -1.116e+01 2.054e+00 3.157e+03 -5.435 5.91e-08 ***
## Time:VisitV2 3.938e-02 1.916e-02 3.152e+03 2.056 0.0399 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Age SexMal MthOnA MthMnU MthDll MthDLM MthDyA Time VistV2
## Age -0.526
## SexMale -0.054 -0.006
## MethOnsetAg -0.570 -0.148 -0.389
## MthMnthsUsd -0.217 -0.386 -0.326 0.575
## MethDollars -0.340 0.120 -0.015 0.172 -0.057
## MthDysLstMn -0.147 -0.114 0.349 -0.022 -0.070 -0.438
## MthDysAbstn -0.505 0.075 0.330 0.054 0.099 0.389 0.054
## Time -0.098 -0.001 -0.006 0.005 0.006 -0.006 -0.003 -0.012
## VisitV2 -0.068 -0.011 0.006 -0.006 0.009 -0.009 0.003 -0.009 0.602
## Time:VistV2 0.065 0.005 0.000 -0.001 -0.005 -0.001 0.005 0.005 -0.695 -0.868
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Age 1696 1696 1 19.06 2.0798 0.16550
## Sex 365 365 1 19.04 0.4474 0.51161
## MethOnsetAge 673 673 1 19.09 0.8249 0.37508
## MethMonthsUsed 86 86 1 18.95 0.1060 0.74828
## MethDollars 9 9 1 19.30 0.0112 0.91663
## MethDaysLastMonth 48 48 1 19.11 0.0589 0.81076
## MethDaysAbstinant 423 423 1 18.97 0.5189 0.48010
## Time 65945 65945 1 3153.62 80.8593 < 2.2e-16 ***
## Visit 24087 24087 1 3156.88 29.5346 5.91e-08 ***
## Time:Visit 3447 3447 1 3152.17 4.2264 0.03988 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
this_lme <- lmer(ImageCompletionTime ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed + OpioidOrHeroinDollars +
OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant + Time * Visit + (1 | id),
data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'opioid',])## Warning: Some predictor variables are on very different scales: consider rescaling
## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'opioid Image Completion Time', show.values = TRUE, type = 'std', order.terms = c(1,2,8, 9, 10, 3,4,6,5,7)))print(summary(this_lme))## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: ImageCompletionTime ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed +
## OpioidOrHeroinDollars + OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant +
## Time * Visit + (1 | id)
## Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==
## "opioid", ]
##
## REML criterion at convergence: 29457.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2662 -0.3552 -0.1292 0.1051 26.5328
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 65.71 8.106
## Residual 598.76 24.470
## Number of obs: 3180, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.465e+01 1.157e+01 1.929e+01 2.996 0.00734 **
## Age 5.175e-02 3.314e-01 1.897e+01 0.156 0.87756
## SexMale 9.412e-01 3.592e+00 1.912e+01 0.262 0.79611
## OpioidOrHeroinOnsetAge 4.239e-01 3.772e-01 1.918e+01 1.124 0.27501
## OpioidOrHeroinMonthsUsed 9.660e-03 2.695e-02 1.906e+01 0.358 0.72398
## OpioidOrHeroinDollars -1.931e-03 1.268e-03 1.916e+01 -1.523 0.14404
## OpioidOrHeroinDaysLastMonth -6.221e-02 2.012e-01 1.895e+01 -0.309 0.76058
## OpioidOrHeroinDaysAbstinant -4.398e-04 1.451e-03 1.899e+01 -0.303 0.76519
## Time -1.401e-01 1.146e-02 3.152e+03 -12.226 < 2e-16 ***
## VisitV2 -1.153e+01 1.734e+00 3.155e+03 -6.648 3.49e-11 ***
## Time:VisitV2 4.732e-02 1.633e-02 3.152e+03 2.898 0.00378 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Age SexMal OpOHOA OpOHMU OpdOHD OOHDLM OpOHDA Time VistV2
## Age -0.752
## SexMale -0.066 0.021
## OpdOrHrnOnA 0.083 -0.616 -0.146
## OpdOrHrnMnU 0.126 -0.396 -0.334 0.407
## OpdOrHrnDll 0.058 0.119 0.212 -0.419 -0.248
## OpdOrHrnDLM -0.631 0.486 0.056 -0.262 -0.393 -0.157
## OpdOrHrnDyA -0.177 -0.045 -0.159 -0.005 0.227 0.017 0.236
## Time -0.084 -0.002 0.003 -0.005 -0.003 0.003 -0.004 0.001
## VisitV2 -0.065 -0.003 0.005 -0.006 -0.006 0.008 -0.004 -0.002 0.603
## Time:VistV2 0.056 0.006 0.002 -0.002 0.001 -0.001 0.004 -0.001 -0.702 -0.863
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Age 15 15 1 18.97 0.0244 0.877562
## Sex 41 41 1 19.12 0.0686 0.796114
## OpioidOrHeroinOnsetAge 756 756 1 19.18 1.2627 0.275008
## OpioidOrHeroinMonthsUsed 77 77 1 19.06 0.1285 0.723983
## OpioidOrHeroinDollars 1389 1389 1 19.16 2.3202 0.144038
## OpioidOrHeroinDaysLastMonth 57 57 1 18.95 0.0956 0.760579
## OpioidOrHeroinDaysAbstinant 55 55 1 18.99 0.0918 0.765194
## Time 121834 121834 1 3151.43 203.4770 < 2.2e-16 ***
## Visit 26463 26463 1 3154.74 44.1971 3.485e-11 ***
## Time:Visit 5029 5029 1 3151.58 8.3997 0.003779 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rc_data1 <- subject_data[subject_data$millisecondsOnPage < 20000 & subject_data$RatingType == 'craving',]
p4 <- ggplot(rc_data1,aes(x=ImageSetCap,y=millisecondsOnPage / 1000, fill = ImageSetCap, colour = ImageSetCap))+
geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim = FALSE)+
geom_point(position = position_jitter(width = .15), size = .25)+
geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = millisecondsOnPage / 1000),outlier.shape = NA,
alpha = 0.3, width = .1, colour = "BLACK") +
ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
colour = FALSE) +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols) +
ggtitle("SecondsOnPage: Craving Ratings") #+ scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))
rc_data2 <- subject_data[subject_data$ImageCompletionTime < 100,]
p5 <- ggplot(rc_data2,aes(x=ImageSetCap,y=ImageCompletionTime, fill = ImageSetCap, colour = ImageSetCap))+
geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim = FALSE)+
geom_point(position = position_jitter(width = .15), size = .25)+
geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = ImageCompletionTime),outlier.shape = NA,
alpha = 0.3, width = .1, colour = "BLACK") +
ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
colour = FALSE) +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols) +
ggtitle("ImageCompletionTime") #+ scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))
png('RTraincloudsExcluded.png', width = 1400, height = 400)
print(annotate_figure(ggarrange(p4, p5, ncol = 1, nrow = 2, common.legend = TRUE, legend = 'right'), top = text_grob('RTRainclouds', color = 'red')))
dev.off()## png
## 2
#remove group summary lines
selection_data <- merged_summary_values[!is.na(merged_summary_values$File),]
#remove opioid image 118 because 109 and 118 are the same
selection_data <- selection_data[(selection_data$Group != 'opioid Slide118.jpeg'),]
#remove images with text on them
selection_data <- selection_data[!(selection_data$Group %in% c('opioid Slide015.jpeg', 'opioid Slide016.jpeg', 'opioid Slide017.jpeg',
'opioid Slide018.jpeg', 'opioid Slide019.jpeg', 'opioid Slide020.jpeg',
'control Slide017.jpeg', 'control Slide018.jpeg', 'control Slide019.jpeg', 'control Slide020.jpeg')),]
#remove light bulb images, since they can be associated with drug use
selection_data <- selection_data[!(selection_data$Group %in% c('opioid Slide57.jpeg', 'opioid Slide87.jpeg')),]
#line below to select data based on hightst/lowest craving
#sorted_data <- selection_data[order(selection_data$allsubjects_craving_mean),]
#instead, let's select based on hue and see if that helps balance things--get images closest to 0.2
try_one_criteria <- function(control_target_hue, control_target_saturation, control_target_value,
meth_target_hue, meth_target_saturation, meth_target_value,
opioid_target_hue, opioid_target_saturation, opioid_target_value){
#sorted_data <- selection_data[sample(nrow(selection_data)),]
selection_data$selection_cost <- NA
selection_data$selection_cost[selection_data$ImageSet == 'control'] <-
(selection_data$hue_mean[selection_data$ImageSet == 'control'] - control_target_hue)^2 +
(selection_data$saturation_mean[selection_data$ImageSet == 'control'] - control_target_saturation)^2 +
(selection_data$value_mean[selection_data$ImageSet == 'control'] - control_target_value)^2
selection_data$selection_cost[selection_data$ImageSet == 'meth'] <-
(selection_data$hue_mean[selection_data$ImageSet == 'meth'] - meth_target_hue)^2 +
(selection_data$saturation_mean[selection_data$ImageSet == 'meth'] - meth_target_saturation)^2 +
(selection_data$value_mean[selection_data$ImageSet == 'meth'] - meth_target_value)^2
selection_data$selection_cost[selection_data$ImageSet == 'opioid'] <-
(selection_data$hue_mean[selection_data$ImageSet == 'opioid'] - opioid_target_hue)^2 +
(selection_data$saturation_mean[selection_data$ImageSet == 'opioid'] - opioid_target_saturation)^2 +
(selection_data$value_mean[selection_data$ImageSet == 'opioid'] - opioid_target_value)^2
(selection_data$hue_mean - 0.3)^2 + (selection_data$saturation_mean - 0.3)^2 + (selection_data$value_mean - 0.3)^2
sorted_data <- selection_data[order(selection_data$selection_cost),]
#label categories 1:6, to make indexind into number_selected easier
sorted_data$category_number <- NA
sorted_data$category_number[sorted_data$Category == 'control_objects'] <- 1
sorted_data$category_number[sorted_data$Category == 'control_objects_hand'] <- 2
sorted_data$category_number[sorted_data$Category == 'control_tool'] <- 3
sorted_data$category_number[sorted_data$Category == 'control_tool_hand_simple'] <- 4
sorted_data$category_number[sorted_data$Category == 'control_tool_hand_complex'] <- 5
sorted_data$category_number[sorted_data$Category == 'control_tool_face'] <- 6
sorted_data$category_number[sorted_data$Category == 'meth'] <- 1
sorted_data$category_number[sorted_data$Category == 'meth_hand'] <- 2
sorted_data$category_number[sorted_data$Category == 'meth_instrument'] <- 3
sorted_data$category_number[sorted_data$Category == 'meth_instrument_hand'] <- 4
sorted_data$category_number[sorted_data$Category == 'meth_injection_hand'] <- 5
sorted_data$category_number[sorted_data$Category == 'meth_face_activities'] <- 6
#not using any meth_injection_instrument images
sorted_data <- sorted_data[sorted_data$Category != 'meth_injection_instrument',]
sorted_data$category_number[sorted_data$Category == 'opioid'] <- 1
sorted_data$category_number[sorted_data$Category == 'opioid_hand'] <- 2
sorted_data$category_number[sorted_data$Category == 'opioid_injection_instrument'] <- 3
sorted_data$category_number[sorted_data$Category == 'opioid_instrument_hand'] <- 4
sorted_data$category_number[sorted_data$Category == 'opioid_injection_hand'] <- 5
sorted_data$category_number[sorted_data$Category == 'opioid_face_activities'] <- 6
#number in each group that have been selected
number_selected <- c(0, 0, 0, 0, 0, 0)
#go forward, marking the first 12 in each category of control images for selection
#will set to TRUE for images selected to be used
#this will be the 12 highest craving drug (opioid and meth separate) and lowest craving neutral images in each category
sorted_data$selected <- FALSE
for (i in 1:nrow(sorted_data)){
#step forwards, taking the least craving-iducing neutral images
if (sorted_data$ImageSet[i] != 'control') {
next
}
this_category <- sorted_data$category_number[i]
if (number_selected[this_category] > 11){
next
}
number_selected[this_category] <- number_selected[this_category] + 1
sorted_data[i, 'selected'] <- TRUE
}
#number in each group that have been selected
number_selected_meth <- c(0, 0, 0, 0, 0, 0)
number_selected_opioid <- c(0, 0, 0, 0, 0, 0)
#step through backwards, selecting most craving-inducing meth/opioid images
#for (i in nrow(sorted_data):1){
for (i in 1:nrow(sorted_data)){
this_category <- sorted_data$category_number[i]
this_group <- sorted_data$ImageSet[i]
if (this_group == 'meth'){
if (number_selected_meth[this_category] > 11){
next
}
number_selected_meth[this_category] <- number_selected_meth[this_category] + 1
sorted_data[i, 'selected'] <- TRUE
}
if (this_group == 'opioid'){
if (number_selected_opioid[this_category] > 11){
next
}
number_selected_opioid[this_category] <- number_selected_opioid[this_category] + 1
sorted_data[i, 'selected'] <- TRUE
}
}
selected_data <- sorted_data[sorted_data$selected,]
set.seed(12345678)
selected_data$in_set <- NA
#will randomize these
in_set <- c(0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2)
for (i in 1:6){
selected_data[selected_data$category_number == i & selected_data$ImageSet == 'control', 'in_set'] <- in_set[sample(12)]
selected_data[selected_data$category_number == i & selected_data$ImageSet == 'meth', 'in_set'] <- in_set[sample(12)]
selected_data[selected_data$category_number == i & selected_data$ImageSet == 'opioid', 'in_set'] <- in_set[sample(12)]
}
library(tableone)
allsets_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean',
'hue_mean', 'saturation_mean', 'value_mean'),
strata = c('ImageSet'), data = selected_data)
print(allsets_tab)
p1 <- as.numeric(substr(print(allsets_tab)[[30]], 2, 6))
p2 <- as.numeric(substr(print(allsets_tab)[[31]], 2, 6))
p3 <- as.numeric(substr(print(allsets_tab)[[32]], 2, 6))
return(selected_data)
}
selected_data <- try_one_criteria(0.4, 0.6, 0.25,
0.3, 0.1, 0.5,
0.2, 0.2, 0.3)## Stratified by ImageSet
## control meth opioid p test
## n 72 72 72
## allsubjects_arousal_mean (mean (sd)) 2.53 (0.25) 4.72 (0.18) 4.65 (0.25) <0.001
## allsubjects_craving_mean (mean (sd)) 11.88 (4.96) 88.54 (3.12) 88.39 (3.68) <0.001
## allsubjects_typicality_mean (mean (sd)) 12.22 (4.55) 89.70 (3.17) 89.54 (3.62) <0.001
## allsubjects_valence_mean (mean (sd)) 4.89 (0.15) 4.42 (0.19) 4.44 (0.25) <0.001
## hue_mean (mean (sd)) 0.24 (0.18) 0.27 (0.14) 0.24 (0.13) 0.275
## saturation_mean (mean (sd)) 0.26 (0.14) 0.30 (0.13) 0.27 (0.13) 0.139
## value_mean (mean (sd)) 0.42 (0.20) 0.39 (0.15) 0.41 (0.15) 0.574
## Stratified by ImageSet
## control meth opioid p test
## n 72 72 72
## allsubjects_arousal_mean (mean (sd)) 2.53 (0.25) 4.72 (0.18) 4.65 (0.25) <0.001
## allsubjects_craving_mean (mean (sd)) 11.88 (4.96) 88.54 (3.12) 88.39 (3.68) <0.001
## allsubjects_typicality_mean (mean (sd)) 12.22 (4.55) 89.70 (3.17) 89.54 (3.62) <0.001
## allsubjects_valence_mean (mean (sd)) 4.89 (0.15) 4.42 (0.19) 4.44 (0.25) <0.001
## hue_mean (mean (sd)) 0.24 (0.18) 0.27 (0.14) 0.24 (0.13) 0.275
## saturation_mean (mean (sd)) 0.26 (0.14) 0.30 (0.13) 0.27 (0.13) 0.139
## value_mean (mean (sd)) 0.42 (0.20) 0.39 (0.15) 0.41 (0.15) 0.574
## Stratified by ImageSet
## control meth opioid p test
## n 72 72 72
## allsubjects_arousal_mean (mean (sd)) 2.53 (0.25) 4.72 (0.18) 4.65 (0.25) <0.001
## allsubjects_craving_mean (mean (sd)) 11.88 (4.96) 88.54 (3.12) 88.39 (3.68) <0.001
## allsubjects_typicality_mean (mean (sd)) 12.22 (4.55) 89.70 (3.17) 89.54 (3.62) <0.001
## allsubjects_valence_mean (mean (sd)) 4.89 (0.15) 4.42 (0.19) 4.44 (0.25) <0.001
## hue_mean (mean (sd)) 0.24 (0.18) 0.27 (0.14) 0.24 (0.13) 0.275
## saturation_mean (mean (sd)) 0.26 (0.14) 0.30 (0.13) 0.27 (0.13) 0.139
## value_mean (mean (sd)) 0.42 (0.20) 0.39 (0.15) 0.41 (0.15) 0.574
## Stratified by ImageSet
## control meth opioid p test
## n 72 72 72
## allsubjects_arousal_mean (mean (sd)) 2.53 (0.25) 4.72 (0.18) 4.65 (0.25) <0.001
## allsubjects_craving_mean (mean (sd)) 11.88 (4.96) 88.54 (3.12) 88.39 (3.68) <0.001
## allsubjects_typicality_mean (mean (sd)) 12.22 (4.55) 89.70 (3.17) 89.54 (3.62) <0.001
## allsubjects_valence_mean (mean (sd)) 4.89 (0.15) 4.42 (0.19) 4.44 (0.25) <0.001
## hue_mean (mean (sd)) 0.24 (0.18) 0.27 (0.14) 0.24 (0.13) 0.275
## saturation_mean (mean (sd)) 0.26 (0.14) 0.30 (0.13) 0.27 (0.13) 0.139
## value_mean (mean (sd)) 0.42 (0.20) 0.39 (0.15) 0.41 (0.15) 0.574
control_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean',
'hue_mean', 'saturation_mean', 'value_mean'),
strata = c('in_set'), data = selected_data[selected_data$ImageSet == 'control',])
control_tab## Stratified by in_set
## 0 1 2 p test
## n 24 24 24
## allsubjects_arousal_mean (mean (sd)) 2.50 (0.26) 2.56 (0.29) 2.54 (0.22) 0.739
## allsubjects_craving_mean (mean (sd)) 11.17 (3.57) 11.61 (3.70) 12.84 (6.91) 0.488
## allsubjects_typicality_mean (mean (sd)) 12.06 (3.37) 11.69 (2.98) 12.92 (6.54) 0.636
## allsubjects_valence_mean (mean (sd)) 4.87 (0.16) 4.90 (0.18) 4.90 (0.13) 0.775
## hue_mean (mean (sd)) 0.28 (0.19) 0.22 (0.18) 0.22 (0.14) 0.394
## saturation_mean (mean (sd)) 0.26 (0.13) 0.25 (0.17) 0.26 (0.11) 0.966
## value_mean (mean (sd)) 0.42 (0.21) 0.39 (0.20) 0.45 (0.20) 0.595
meth_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean',
'hue_mean', 'saturation_mean', 'value_mean'),
strata = c('in_set'), data = selected_data[selected_data$ImageSet == 'meth',])
meth_tab## Stratified by in_set
## 0 1 2 p test
## n 24 24 24
## allsubjects_arousal_mean (mean (sd)) 4.71 (0.22) 4.74 (0.17) 4.70 (0.17) 0.715
## allsubjects_craving_mean (mean (sd)) 87.35 (2.71) 89.26 (3.60) 89.00 (2.73) 0.070
## allsubjects_typicality_mean (mean (sd)) 88.34 (2.91) 90.30 (3.33) 90.46 (2.93) 0.034
## allsubjects_valence_mean (mean (sd)) 4.42 (0.25) 4.40 (0.15) 4.44 (0.17) 0.817
## hue_mean (mean (sd)) 0.25 (0.12) 0.30 (0.18) 0.27 (0.12) 0.395
## saturation_mean (mean (sd)) 0.29 (0.14) 0.30 (0.13) 0.32 (0.12) 0.683
## value_mean (mean (sd)) 0.36 (0.11) 0.43 (0.17) 0.38 (0.15) 0.204
opioid_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean',
'hue_mean', 'saturation_mean', 'value_mean'),
strata = c('in_set'), data = selected_data[selected_data$ImageSet == 'opioid',])
opioid_tab## Stratified by in_set
## 0 1 2 p test
## n 24 24 24
## allsubjects_arousal_mean (mean (sd)) 4.69 (0.23) 4.65 (0.25) 4.62 (0.26) 0.624
## allsubjects_craving_mean (mean (sd)) 88.86 (3.29) 88.72 (4.17) 87.60 (3.56) 0.439
## allsubjects_typicality_mean (mean (sd)) 90.07 (2.99) 89.32 (4.52) 89.23 (3.26) 0.684
## allsubjects_valence_mean (mean (sd)) 4.41 (0.23) 4.48 (0.27) 4.44 (0.27) 0.575
## hue_mean (mean (sd)) 0.23 (0.15) 0.24 (0.10) 0.24 (0.12) 0.865
## saturation_mean (mean (sd)) 0.25 (0.12) 0.28 (0.13) 0.30 (0.15) 0.556
## value_mean (mean (sd)) 0.39 (0.13) 0.43 (0.17) 0.41 (0.16) 0.660
set0_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean',
'hue_mean', 'saturation_mean', 'value_mean'),
strata = c('ImageSet'), data = selected_data[selected_data$in_set == 0,])
print(set0_tab)## Stratified by ImageSet
## control meth opioid p test
## n 24 24 24
## allsubjects_arousal_mean (mean (sd)) 2.50 (0.26) 4.71 (0.22) 4.69 (0.23) <0.001
## allsubjects_craving_mean (mean (sd)) 11.17 (3.57) 87.35 (2.71) 88.86 (3.29) <0.001
## allsubjects_typicality_mean (mean (sd)) 12.06 (3.37) 88.34 (2.91) 90.07 (2.99) <0.001
## allsubjects_valence_mean (mean (sd)) 4.87 (0.16) 4.42 (0.25) 4.41 (0.23) <0.001
## hue_mean (mean (sd)) 0.28 (0.19) 0.25 (0.12) 0.23 (0.15) 0.515
## saturation_mean (mean (sd)) 0.26 (0.13) 0.29 (0.14) 0.25 (0.12) 0.604
## value_mean (mean (sd)) 0.42 (0.21) 0.36 (0.11) 0.39 (0.13) 0.451
set1_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean',
'hue_mean', 'saturation_mean', 'value_mean'),
strata = c('ImageSet'), data = selected_data[selected_data$in_set == 1,])
print(set1_tab)## Stratified by ImageSet
## control meth opioid p test
## n 24 24 24
## allsubjects_arousal_mean (mean (sd)) 2.56 (0.29) 4.74 (0.17) 4.65 (0.25) <0.001
## allsubjects_craving_mean (mean (sd)) 11.61 (3.70) 89.26 (3.60) 88.72 (4.17) <0.001
## allsubjects_typicality_mean (mean (sd)) 11.69 (2.98) 90.30 (3.33) 89.32 (4.52) <0.001
## allsubjects_valence_mean (mean (sd)) 4.90 (0.18) 4.40 (0.15) 4.48 (0.27) <0.001
## hue_mean (mean (sd)) 0.22 (0.18) 0.30 (0.18) 0.24 (0.10) 0.186
## saturation_mean (mean (sd)) 0.25 (0.17) 0.30 (0.13) 0.28 (0.13) 0.567
## value_mean (mean (sd)) 0.39 (0.20) 0.43 (0.17) 0.43 (0.17) 0.714
set2_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean',
'hue_mean', 'saturation_mean', 'value_mean'),
strata = c('ImageSet'), data = selected_data[selected_data$in_set == 2,])
print(set2_tab)## Stratified by ImageSet
## control meth opioid p test
## n 24 24 24
## allsubjects_arousal_mean (mean (sd)) 2.54 (0.22) 4.70 (0.17) 4.62 (0.26) <0.001
## allsubjects_craving_mean (mean (sd)) 12.84 (6.91) 89.00 (2.73) 87.60 (3.56) <0.001
## allsubjects_typicality_mean (mean (sd)) 12.92 (6.54) 90.46 (2.93) 89.23 (3.26) <0.001
## allsubjects_valence_mean (mean (sd)) 4.90 (0.13) 4.44 (0.17) 4.44 (0.27) <0.001
## hue_mean (mean (sd)) 0.22 (0.14) 0.27 (0.12) 0.24 (0.12) 0.402
## saturation_mean (mean (sd)) 0.26 (0.11) 0.32 (0.12) 0.30 (0.15) 0.296
## value_mean (mean (sd)) 0.45 (0.20) 0.38 (0.15) 0.41 (0.16) 0.333
write.csv(selected_data, 'DCR_selected_split-MatchedHSV.csv', row.names=FALSE, quote = FALSE)scatter_data <- subject_data[subject_data$RatingType == c('craving') & (subject_data$millisecondsOnPage < 20000) & subject_data$ImageSet %in% c('meth', 'opioid'), ]
ggplot(scatter_data, aes(x = Value, y = millisecondsOnPage / 1000, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm') +
xlab('Craving') + ylab('Seconds') +
ggtitle(paste0('Craving vs Craving RT')) +
scale_color_manual(values = cols)scatter_data <- subject_data[subject_data$RatingType == c('craving') & (subject_data$ImageCompletionTime < 100) & subject_data$ImageSet %in% c('meth', 'opioid'), ]
ggplot(scatter_data, aes(x = Value, y = ImageCompletionTime, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm') +
xlab('Craving') + ylab('Seconds') +
ggtitle(paste0('Craving vs Time On Image')) +
scale_color_manual(values = cols)